Analyze.rmd

---

Description:

  This script takes inputs in the form of counts matrices and dataframes   formatted for analysis by the CountstoAnalysis scripts for the 001 and   002 datasets (with all calculations and formatting already completed) –   See Inputs for all inputs. The script then conducts all summary   calculations, statistical tests, and produces all visuals and output tables   found in the main text and supplement (unless otherwise noted below – see   Outputs). Code blocks appear in more-or-less the same order that data are   presented in the text; brief descriptions head each box and more detailed   code can be found in-line. This script outputs a set of figures and tables   that are included, some with modification and stylization, in the main text   and supplement.

---

Inputs:

1. poc1.Rdata – CountstoAnalysis formatted dataframe for POC fitness assay single-barcode controls.
2. poc2.Rdata – CountstoAnalysis formatted dataframe for POC fitness assay dual-barcode controls.
3. poc92.Rdata – CountstoAnalysis formatted dataframe for POC fitness assay fitness assay samples.
4. 001_counts.expected.Rdata – FastqtoCounts formatted matrix for POC fitness assay library expected (by design) counts.
5. 001_counts.unexpected.Rdata – FastqtoCounts formatted matrix for POC fitness assay library unexpected (barcode cross-contamination) counts.
6. 002_CC_counts.expected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment expected (by design) counts – consensus matrix constructed from two sequenced libraries.
7. 002_CC_counts.unexpected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment unexpected (barcode cross-contamination) counts – consensus matrix from two sequenced libraries.
8. 002_counts.expected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment expected (by design) counts – original sequenced library.
9. 002_counts.unexpected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment unexpected (barcode cross-contamination) counts – original sequenced library.
10. 002_Reruns_counts.expected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment expected (by design) counts – resequenced samples library.
11. 002_Reruns_counts.unexpected.Rdata – FastqtoCounts formatted matrix 250-generation evolution experiment unexpected (barcode cross-contamination) counts – resequenced samples library.
12. myfa.Rdata – CountstoAnalysis formatted dataframe for 250-generation experimental evolution project fitness assays.
13. myevo.Rdata – CountstoAnalysis formatted dataframe for 250-generation experimental evolution project longitudinal samples.

---

Outputs:

1. 000_Figure_2.pdf – Manuscript Figure 2 output as .pdf file.
   
2. 000_Figure_3.pdf – Manuscript Figure 3 output as .pdf file.
3. 000_Figure_4.pdf – Manuscript Figure 4 output as .pdf file.
4. 000_Figure_S1_A.pdf – Manuscript Figure S1 panel a. output as .pdf file.
5. 000_Figure_S1_B.pdf – Manuscript Figure S1 panel b. output as .pdf file.
6. 000_Figure_S2.pdf – Manuscript Figure S2 output as .pdf file.
7. 000_Figure_S3.pdf – Manuscript Figure S3 output as .pdf file.
8. 000_Figure_S4.pdf – Manuscript Figure S4 output as .pdf file.
9. 000_Figure_S5.pdf – Manuscript Figure S5 output as .pdf file.
10. 000_Figure_S6.pdf – Manuscript Figure S6 output as .pdf file.
11. 000_Figure_S7.pdf – Manuscript Figure S7 output as .pdf file.
12. 000_Figure_S8.pdf – Manuscript Figure S8 output as .pdf file.
13. 000_Figure_S9.pdf – Manuscript Figure S9 output as .pdf file.
14. 000_Figure_S10.pdf – Manuscript Figure S10 output as .pdf file.
15. 000_Table_S4.html – Raw Manuscript Figure S4 output as .html file. File later stylized in excel and P-values converted to fdr for publication.
16. 000_Table_S5.csv – Raw Manuscript Figure S5 output as .csv file. File later stylized in excel, additional columns added.
17. 000_Table_S6.html – Raw Manuscript Figure S6 output as .html file. File later stylized in excel.
18. 000_Table_S7.html – Raw Manuscript Figure S7 output as .html file. File later stylized in excel.
19. 000_Table_S8.html – Raw Manuscript Figure S8 output as .html file. File later stylized in excel and P-values converted to 1-sided (higher) and fdr for publication.
20. 000_Table_S9.html – Raw Manuscript Figure S9 output as .html file. File later stylized in excel.
21. 000_Table_S10.html – Raw Manuscript Figure S10 output as .html file. File later stylized in excel.
22. 000_Table_S11.html – Raw Manuscript Figure S11 output as .html file. File later stylized in excel.
23. 000_Table_S12.html – Raw Manuscript Figure S12 output as .html file. File later stylized in excel.
24. 000_Table_S13.html – Raw Manuscript Figure S13 output as .html file. File later stylized in excel.
25. 000_Table_S14.html – Raw Manuscript Figure S14 output as .html file. File later stylized in excel.
26. 000_Table_S15.html – Raw Manuscript Figure S15 output as .html file. File later stylized in excel.

• Table_S4 through Table_S16 are all stylized and formatted in Excel and converted to .xlsx files from .html (if originally .html formatted) for publication. ** Table_S4 is output with raw p-values, but values are converted to fdr for publication. Table_S8 is output with raw 2-sided test p-values, but values are converted to 1-sided (higher) and fdr for publication. *** Figure_1, Figure_5 are not output by this script. They are manually created in Keynote and Excel, respectively. **** Table_S1, Table_S2, and Table_S3 are not output by this script. They are manually created in Excel.

---

Begin Code Blocks:

Prepare the workspace: clean, set options, load packages, load and final formatting on inputs.

rm(list = ls()) # clear workspace
knitr::opts_chunk$set(tidy = TRUE, collapse = TRUE) # set global knitr options.
# load required packages
require(ggplot2); # used for plotting
require(pwr); # used for power analysis
require(lme4); # used for mixed effects models
require(lmerTest); # used to interpret mixed effects models
require(sjPlot); # used to output tables
require(weights); # used for weighted t.test
require(plyr); # used for dataframe manipulation

options(scipen = 999) # turn off scientific notation

getwd()
## [1] "/Users/boterolab1/Box Sync/CB_VF_Shared/Wet_Lab/Projects/MP_Clean/Methods_Manuscript_V6/File_S3/3_Analysis"

dir_in_pocpro <- "2_CountstoAnalysis/001_Prepare_Datasets"  # input directory for Proof of concept processed dataframes
setwd("..")
setwd(dir_in_pocpro)  # set to above directory path
load(file = "poc1.Rdata")  # load the 1bc control processed dataframe
load(file = "poc2.Rdata")  # load the 2 bc control processed dataframe
load(file = "poc92.Rdata")  # load the 92 bc proof of concept fitness assays processed dataframe

dir_in_pocraw <- "001_Raw_Data_Metadata"  # input directory for Proof of concept raw dataframes
setwd("..")
setwd(dir_in_pocraw)  # set to above directory path
load(file = "001_counts.expected.Rdata")
poc.ce <- (counts.expected + 1)[c(1:92), ]
rm(counts.expected)  # expected counts
load(file = "001_counts.unexpected.Rdata")
poc.cc <- (counts.unexpected + 1)[c(1:92), ]
rm(counts.unexpected)  # cross contamination counts. 

dir_in_faevoraw <- "002_Raw_Data_Metadata"  # input directory for FA and Evo raw dataframes
setwd("..")
setwd(dir_in_faevoraw)  # set to above directory path
load(file = "002_CC_counts.expected.Rdata")
faevo.ce <- (ects + 0)
rm(ects)  # expected counts {combined}
load(file = "002_CC_counts.unexpected.Rdata")
faevo.cc <- (ccts + 0)
rm(ccts)  # cross contamination counts.{combined}
load(file = "002_counts.expected.Rdata")
faevo.1.ce <- (counts.expected + 1)[c(1:78), ]
rm(counts.expected)  # expected counts {first run samples}
load(file = "002_counts.unexpected.Rdata")
faevo.1.cc <- (counts.unexpected + 1)[c(1:78), ]
rm(counts.unexpected)  # cross contamination counts. {first run samples}
load(file = "002_Reruns_counts.expected.Rdata")
faevo.2.ce <- (counts.expected + 1)[c(1:78), ]
rm(counts.expected)  # expected counts {rerun samples}
load(file = "002_Reruns_counts.unexpected.Rdata")
faevo.2.cc <- (counts.unexpected + 1)[c(1:78), ]
rm(counts.unexpected)  # cross contamination counts. {rerun samples}

dir_in_faevopro <- "002_Prepare_Datasets"  # input directory for FA and Evo processed dataframes
setwd("..")
setwd(dir_in_faevopro)  # set to above directory path
load(file = "myfa.Rdata")  # fitness assay
load(file = "myevo.Rdata")  # evolution
rm(dir_in_faevopro, dir_in_faevoraw, dir_in_pocpro, dir_in_pocraw)

Create helper functions for linear model residual plotting and data transformation.

# expand_residuals creates a modified histogram of residuals. Each entry in
# the model output is plotted a number of times equal to its reads.
expand_residuals <- function(model, reads) {
    myfr <- NA
    for (i in 1:length(residuals(model))) {
        times <- (((reads[i] - min(reads, na.rm = T))/(max(reads) - min(reads))) * 
            1000) + 1
        myfr <- c(myfr, rep(residuals(model)[i], times = times))
    }
    return(myfr)
}

# log transformation for modeling.
log_transform <- function(myvar) {
    log(((myvar - min(myvar))/(max(myvar) - min(myvar))) + 1)
    return(myvar)
}

Set color palette options for all plots

setpalette <- RColorBrewer::brewer.pal(6, "Dark2")  # current preferred palette
# setpalette <-
# c('#993404','#d95f0e','#fe9929','#41b6c4','#2c7fb8','#253494') # option 2
# setpalette <- c('#8c510a','#d8b365','#f6e8c3',
# '#c7eae5','#5ab4ac','#01665e') # option 3

Analysis: Assess POC fitness assay barcode fitness equivalence.

my <- poc92[order(poc92$bcid), ]

# model & diagnostic plots
# ----------------------------------------------------
my$rw_corrected <- my$rw - weighted.mean(my$rw, my$r)
mymod <- lm(my$rw_corrected ~ my$bcid + 0, weights = my$r)
summary(mymod)
## 
## Call:
## lm(formula = my$rw_corrected ~ my$bcid + 0, weights = my$r)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7350 -0.5858 -0.0109  0.5304  5.3405 
## 
## Coefficients:
##                 Estimate  Std. Error t value             Pr(>|t|)    
## my$bcidd1A10  0.00132555  0.00371437   0.357             0.721280    
## my$bcidd1A11 -0.00762370  0.00466055  -1.636             0.102267    
## my$bcidd1A2   0.00811450  0.00396019   2.049             0.040778 *  
## my$bcidd1A3   0.00494242  0.00439553   1.124             0.261165    
## my$bcidd1A4   0.00795573  0.00448555   1.774             0.076495 .  
## my$bcidd1A5   0.00103320  0.00520070   0.199             0.842573    
## my$bcidd1A6  -0.00085773  0.00357809  -0.240             0.810610    
## my$bcidd1A7   0.00344369  0.00512097   0.672             0.501476    
## my$bcidd1A8   0.00583557  0.00414724   1.407             0.159778    
## my$bcidd1A9   0.00724382  0.00374750   1.933             0.053584 .  
## my$bcidd1B1  -0.00527477  0.00502129  -1.050             0.293807    
## my$bcidd1B10 -0.00821202  0.00615734  -1.334             0.182674    
## my$bcidd1B11 -0.01330702  0.00511428  -2.602             0.009437 ** 
## my$bcidd1B12 -0.00741678  0.00553232  -1.341             0.180413    
## my$bcidd1B2   0.00289719  0.00431315   0.672             0.501957    
## my$bcidd1B3   0.00354456  0.00434513   0.816             0.414878    
## my$bcidd1B4   0.00951091  0.00421461   2.257             0.024293 *  
## my$bcidd1B5   0.00961576  0.00421606   2.281             0.022819 *  
## my$bcidd1B6   0.00134667  0.00501745   0.268             0.788462    
## my$bcidd1B7  -0.00398474  0.00435334  -0.915             0.360287    
## my$bcidd1B8   0.00386754  0.00435584   0.888             0.374856    
## my$bcidd1B9  -0.00197406  0.00517408  -0.382             0.702909    
## my$bcidd1C1  -0.00130215  0.00523917  -0.249             0.803778    
## my$bcidd1C10 -0.00200878  0.00612990  -0.328             0.743221    
## my$bcidd1C11  0.00332553  0.00513862   0.647             0.517707    
## my$bcidd1C12 -0.00495638  0.00525931  -0.942             0.346266    
## my$bcidd1C2   0.00310014  0.00440966   0.703             0.482235    
## my$bcidd1C3   0.00170870  0.00544151   0.314             0.753591    
## my$bcidd1C4   0.00064326  0.00445091   0.145             0.885124    
## my$bcidd1C5  -0.00293377  0.00470909  -0.623             0.533457    
## my$bcidd1C6  -0.00108375  0.00401081  -0.270             0.787069    
## my$bcidd1C7  -0.00224296  0.00857907  -0.261             0.793814    
## my$bcidd1C8  -0.00815331  0.00535937  -1.521             0.128566    
## my$bcidd1C9   0.01983075  0.00598788   3.312             0.000968 ***
## my$bcidd1D1   0.00371456  0.00415556   0.894             0.371651    
## my$bcidd1D10  0.01036458  0.00405675   2.555             0.010802 *  
## my$bcidd1D11 -0.00004981  0.00380942  -0.013             0.989570    
## my$bcidd1D12  0.00301900  0.00386910   0.780             0.435449    
## my$bcidd1D2  -0.01098378  0.02236414  -0.491             0.623464    
## my$bcidd1D3  -0.00440055  0.00562205  -0.783             0.434011    
## my$bcidd1D4  -0.00693858  0.00737022  -0.941             0.346760    
## my$bcidd1D5  -0.00326940  0.00563709  -0.580             0.562087    
## my$bcidd1D6  -0.00844038  0.00554061  -1.523             0.128053    
## my$bcidd1D7  -0.00168313  0.00462568  -0.364             0.716052    
## my$bcidd1D8   0.00180750  0.00397535   0.455             0.649462    
## my$bcidd1D9   0.00014080  0.00545635   0.026             0.979419    
## my$bcidd1E1  -0.01270595  0.00523534  -2.427             0.015441 *  
## my$bcidd1E10 -0.00147548  0.00433537  -0.340             0.733691    
## my$bcidd1E11 -0.00708010  0.00406068  -1.744             0.081609 .  
## my$bcidd1E12  0.01157108  0.00442433   2.615             0.009078 ** 
## my$bcidd1E2   0.00500682  0.00373187   1.342             0.180085    
## my$bcidd1E3   0.00762603  0.00554298   1.376             0.169259    
## my$bcidd1E4   0.00591987  0.00384425   1.540             0.123964    
## my$bcidd1E5  -0.00657763  0.00498803  -1.319             0.187644    
## my$bcidd1E6  -0.00139332  0.00372470  -0.374             0.708444    
## my$bcidd1E7   0.00852322  0.00396829   2.148             0.032020 *  
## my$bcidd1E8  -0.00369015  0.00443221  -0.833             0.405327    
## my$bcidd1E9   0.00639603  0.00436339   1.466             0.143076    
## my$bcidd1F1  -0.00012498  0.00381935  -0.033             0.973904    
## my$bcidd1F10  0.00393912  0.00423405   0.930             0.352467    
## my$bcidd1F11  0.00670256  0.00469224   1.428             0.153547    
## my$bcidd1F12  0.00778245  0.00385357   2.020             0.043756 *  
## my$bcidd1F2   0.00582870  0.00435141   1.339             0.180780    
## my$bcidd1F3   0.00651265  0.00596554   1.092             0.275281    
## my$bcidd1F4  -0.02625366  0.00923796  -2.842             0.004596 ** 
## my$bcidd1F5  -0.04834088  0.00429859 -11.246 < 0.0000000000000002 ***
## my$bcidd1F6  -0.00012930  0.00381509  -0.034             0.972972    
## my$bcidd1F7  -0.01483489  0.00492987  -3.009             0.002700 ** 
## my$bcidd1F8  -0.00220117  0.00605837  -0.363             0.716455    
## my$bcidd1F9   0.00029767  0.01166664   0.026             0.979651    
## my$bcidd1G1   0.01006985  0.00654635   1.538             0.124376    
## my$bcidd1G10  0.00690433  0.00358418   1.926             0.054408 .  
## my$bcidd1G11  0.00299309  0.00557059   0.537             0.591205    
## my$bcidd1G12  0.00147423  0.00344561   0.428             0.668867    
## my$bcidd1G2   0.00615615  0.00569582   1.081             0.280096    
## my$bcidd1G3   0.00097634  0.00390830   0.250             0.802795    
## my$bcidd1G4  -0.00335673  0.00478773  -0.701             0.483433    
## my$bcidd1G5  -0.01510445  0.00537164  -2.812             0.005043 ** 
## my$bcidd1G6  -0.00650477  0.00473025  -1.375             0.169463    
## my$bcidd1G7  -0.00911754  0.00430618  -2.117             0.034534 *  
## my$bcidd1G8  -0.00467010  0.00489124  -0.955             0.339967    
## my$bcidd1G9  -0.00345997  0.00470652  -0.735             0.462463    
## my$bcidd1H11  0.00188567  0.00387909   0.486             0.627018    
## my$bcidd1H2  -0.00688077  0.00535513  -1.285             0.199193    
## my$bcidd1H3  -0.00074760  0.00640698  -0.117             0.907138    
## my$bcidd1H4  -0.03016766  0.00925002  -3.261             0.001155 ** 
## my$bcidd1H5  -0.00289320  0.00438300  -0.660             0.509377    
## my$bcidd1H6   0.00117090  0.00365259   0.321             0.748621    
## my$bcidd1H7  -0.00459377  0.00852914  -0.539             0.590312    
## my$bcidd1H8  -0.00449698  0.00614172  -0.732             0.464255    
## my$bcidd1H9  -0.00798745  0.00584794  -1.366             0.172359    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.127 on 819 degrees of freedom
## Multiple R-squared:  0.2653, Adjusted R-squared:  0.1836 
## F-statistic:  3.25 on 91 and 819 DF,  p-value: < 0.00000000000000022
hist(residuals(mymod), breaks = 50)  # residuals look good, check the other diagnostic plots too...

hist(expand_residuals(mymod, my$r), breaks = 50)  # this more accurately represents the distribution of residuals in the model output.

plot(mymod)  # view the diagnostic plots

paste0("Number of BCs with fitness different than pop mean = ", sum(p.adjust(coef(summary(mymod))[, 
    4], method = "fdr") < 0.05))  # number of barcodes with fit different than pop weighted mean. 
## [1] "Number of BCs with fitness different than pop mean = 3"
paste0("RsMSE = ", sqrt(mean(residuals(mymod)^2)))  # report rmse (requires unscalled lm model)
## [1] "RsMSE = 0.0175768866199638"
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S4.html")

	my$rw corrected
Predictors	Estimates	CI	Statistic	p
d 1 A 10	0.00133	-0.00595 – 0.00861	0.35687	0.72127977
d 1 A 11	-0.00762	-0.01676 – 0.00151	-1.63579	0.10226678
d 1 A 2	0.00811	0.00035 – 0.01588	2.04902	0.04077841
d 1 A 3	0.00494	-0.00367 – 0.01356	1.12442	0.26116464
d 1 A 4	0.00796	-0.00084 – 0.01675	1.77364	0.07649506
d 1 A 5	0.00103	-0.00916 – 0.01123	0.19867	0.84257323
d 1 A 6	-0.00086	-0.00787 – 0.00616	-0.23972	0.81061012
d 1 A 7	0.00344	-0.00659 – 0.01348	0.67247	0.50147576
d 1 A 8	0.00584	-0.00229 – 0.01396	1.40710	0.15977841
d 1 A 9	0.00724	-0.00010 – 0.01459	1.93297	0.05358384
d 1 B 1	-0.00527	-0.01512 – 0.00457	-1.05048	0.29380661
d 1 B 10	-0.00821	-0.02028 – 0.00386	-1.33370	0.18267430
d 1 B 11	-0.01331	-0.02333 – -0.00328	-2.60193	0.00943744
d 1 B 12	-0.00742	-0.01826 – 0.00343	-1.34063	0.18041306
d 1 B 2	0.00290	-0.00556 – 0.01135	0.67171	0.50195688
d 1 B 3	0.00354	-0.00497 – 0.01206	0.81575	0.41487762
d 1 B 4	0.00951	0.00125 – 0.01777	2.25665	0.02429251
d 1 B 5	0.00962	0.00135 – 0.01788	2.28075	0.02281947
d 1 B 6	0.00135	-0.00849 – 0.01118	0.26840	0.78846186
d 1 B 7	-0.00398	-0.01252 – 0.00455	-0.91533	0.36028749
d 1 B 8	0.00387	-0.00467 – 0.01240	0.88790	0.37485596
d 1 B 9	-0.00197	-0.01212 – 0.00817	-0.38153	0.70290904
d 1 C 1	-0.00130	-0.01157 – 0.00897	-0.24854	0.80377774
d 1 C 10	-0.00201	-0.01402 – 0.01001	-0.32770	0.74322104
d 1 C 11	0.00333	-0.00675 – 0.01340	0.64716	0.51770702
d 1 C 12	-0.00496	-0.01526 – 0.00535	-0.94240	0.34626553
d 1 C 2	0.00310	-0.00554 – 0.01174	0.70303	0.48223454
d 1 C 3	0.00171	-0.00896 – 0.01237	0.31401	0.75359137
d 1 C 4	0.00064	-0.00808 – 0.00937	0.14452	0.88512380
d 1 C 5	-0.00293	-0.01216 – 0.00630	-0.62300	0.53345671
d 1 C 6	-0.00108	-0.00894 – 0.00678	-0.27021	0.78706863
d 1 C 7	-0.00224	-0.01906 – 0.01457	-0.26145	0.79381434
d 1 C 8	-0.00815	-0.01866 – 0.00235	-1.52132	0.12856573
d 1 C 9	0.01983	0.00809 – 0.03157	3.31181	0.00096754
d 1 D 1	0.00371	-0.00443 – 0.01186	0.89388	0.37165067
d 1 D 10	0.01036	0.00241 – 0.01832	2.55490	0.01080165
d 1 D 11	-0.00005	-0.00752 – 0.00742	-0.01308	0.98957011
d 1 D 12	0.00302	-0.00456 – 0.01060	0.78029	0.43544853
d 1 D 2	-0.01098	-0.05482 – 0.03285	-0.49113	0.62346368
d 1 D 3	-0.00440	-0.01542 – 0.00662	-0.78273	0.43401108
d 1 D 4	-0.00694	-0.02138 – 0.00751	-0.94144	0.34675951
d 1 D 5	-0.00327	-0.01432 – 0.00778	-0.57998	0.56208732
d 1 D 6	-0.00844	-0.01930 – 0.00242	-1.52337	0.12805285
d 1 D 7	-0.00168	-0.01075 – 0.00738	-0.36387	0.71605183
d 1 D 8	0.00181	-0.00598 – 0.00960	0.45468	0.64946246
d 1 D 9	0.00014	-0.01055 – 0.01084	0.02581	0.97941897
d 1 E 1	-0.01271	-0.02297 – -0.00244	-2.42696	0.01544074
d 1 E 10	-0.00148	-0.00997 – 0.00702	-0.34034	0.73369072
d 1 E 11	-0.00708	-0.01504 – 0.00088	-1.74357	0.08160859
d 1 E 12	0.01157	0.00290 – 0.02024	2.61533	0.00907808
d 1 E 2	0.00501	-0.00231 – 0.01232	1.34164	0.18008470
d 1 E 3	0.00763	-0.00324 – 0.01849	1.37580	0.16925925
d 1 E 4	0.00592	-0.00161 – 0.01345	1.53993	0.12396449
d 1 E 5	-0.00658	-0.01635 – 0.00320	-1.31868	0.18764413
d 1 E 6	-0.00139	-0.00869 – 0.00591	-0.37408	0.70844388
d 1 E 7	0.00852	0.00075 – 0.01630	2.14783	0.03201970
d 1 E 8	-0.00369	-0.01238 – 0.00500	-0.83258	0.40532715
d 1 E 9	0.00640	-0.00216 – 0.01495	1.46584	0.14307587
d 1 F 1	-0.00012	-0.00761 – 0.00736	-0.03272	0.97390392
d 1 F 10	0.00394	-0.00436 – 0.01224	0.93034	0.35246715
d 1 F 11	0.00670	-0.00249 – 0.01590	1.42844	0.15354746
d 1 F 12	0.00778	0.00023 – 0.01534	2.01954	0.04375588
d 1 F 2	0.00583	-0.00270 – 0.01436	1.33950	0.18078013
d 1 F 3	0.00651	-0.00518 – 0.01820	1.09171	0.27528055
d 1 F 4	-0.02625	-0.04436 – -0.00815	-2.84193	0.00459550
d 1 F 5	-0.04834	-0.05677 – -0.03992	-11.24575	<0.001
d 1 F 6	-0.00013	-0.00761 – 0.00735	-0.03389	0.97297191
d 1 F 7	-0.01483	-0.02450 – -0.00517	-3.00919	0.00269975
d 1 F 8	-0.00220	-0.01408 – 0.00967	-0.36333	0.71645463
d 1 F 9	0.00030	-0.02257 – 0.02316	0.02551	0.97965060
d 1 G 1	0.01007	-0.00276 – 0.02290	1.53824	0.12437600
d 1 G 10	0.00690	-0.00012 – 0.01393	1.92634	0.05440800
d 1 G 11	0.00299	-0.00793 – 0.01391	0.53730	0.59120517
d 1 G 12	0.00147	-0.00528 – 0.00823	0.42786	0.66886664
d 1 G 2	0.00616	-0.00501 – 0.01732	1.08082	0.28009589
d 1 G 3	0.00098	-0.00668 – 0.00864	0.24981	0.80279516
d 1 G 4	-0.00336	-0.01274 – 0.00603	-0.70111	0.48343252
d 1 G 5	-0.01510	-0.02563 – -0.00458	-2.81189	0.00504284
d 1 G 6	-0.00650	-0.01578 – 0.00277	-1.37514	0.16946284
d 1 G 7	-0.00912	-0.01756 – -0.00068	-2.11732	0.03453391
d 1 G 8	-0.00467	-0.01426 – 0.00492	-0.95479	0.33996652
d 1 G 9	-0.00346	-0.01268 – 0.00576	-0.73514	0.46246280
d 1 H 11	0.00189	-0.00572 – 0.00949	0.48611	0.62701817
d 1 H 2	-0.00688	-0.01738 – 0.00362	-1.28489	0.19919285
d 1 H 3	-0.00075	-0.01331 – 0.01181	-0.11669	0.90713761
d 1 H 4	-0.03017	-0.04830 – -0.01204	-3.26136	0.00115463
d 1 H 5	-0.00289	-0.01148 – 0.00570	-0.66010	0.50937740
d 1 H 6	0.00117	-0.00599 – 0.00833	0.32057	0.74862140
d 1 H 7	-0.00459	-0.02131 – 0.01212	-0.53860	0.59031160
d 1 H 8	-0.00450	-0.01653 – 0.00754	-0.73220	0.46425473
d 1 H 9	-0.00799	-0.01945 – 0.00347	-1.36586	0.17235858
Observations	910
R2 / R2 adjusted	0.265 / 0.184

# -----------------------------------------------------------------------------

# create temporary plotting dataframe
# -----------------------------------------
myres <- as.data.frame(coef(summary(mymod)), stringsAsFactors = F)  # create a myres dataframe to aid in plotting. 
myres$p.adjust <- p.adjust(myres$`Pr(>|t|)`, method = "fdr")
myres$bcid <- gsub("^.*?d", "", rownames(myres))  # fix barcode id names
myres$mrw <- ddply(my, ~bcid, function(my) weighted.mean(my$rw, my$r))[, 2]
myres$rmrw <- ddply(my, ~bcid, function(my) mean(my$r))[, 2]
myres$mcmrw <- myres$mrw - weighted.mean(myres$mrw, myres$rmrw)
myres$color <- "fdr > 0.05"
myres$color[myres$p.adjust <= 0.05] <- "fdr <= 0.05"
myres$color <- as.factor(myres$color)
myres$color <- relevel(myres$color, "fdr > 0.05")
# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- myres$mcmrw
myr <- myres$rmrw
myc <- myres$color
mybinwidth <- 0.0025

mytitle <- NULL
myylab <- "Count"
myxlab <- "Fitness difference from population mean"

myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette

mytextsize <- 8
myheight <- 3.3
mywidth <- 3.3

p <- ggplot(myres, aes(x = myx, color = myc, fill = myc)) + geom_histogram(alpha = myfillapha, 
    position = "identity", binwidth = mybinwidth) + geom_rug(alpha = myptalpha) + 
    xlim(c(-0.06, 0.06)) + ylim(c(NA, 15)) + geom_vline(xintercept = 0, linetype = "dashed", 
    color = myhighlight) + scale_color_manual(values = mypalette, guide = FALSE) + 
    scale_fill_manual(values = mypalette) + theme_classic() + theme(legend.position = c(0.8, 
    0.8)) + labs(fill = "Significance") + ggtitle(mytitle) + ylab(myylab) + 
    xlab(myxlab) + geom_hline(yintercept = -0.01, colour = "white", size = 1.01) + 
    theme(axis.title.y = element_text(size = mytextsize, face = "bold", margin = margin(0, 
        0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) + 
    theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p
## Warning: Removed 4 rows containing missing values (geom_bar).

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_2.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 4 rows containing missing values (geom_bar).
dev.off()  # one column wide, equal height
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(p, myres, mymod, my)

Analysis: Calculate sequenced library reads summary values

paste0("matching counts (match f and r primer barcode, and moby bc): ", sum(c(poc.ce, 
    poc.cc, faevo.1.ce, faevo.1.cc, faevo.2.ce, faevo.2.cc), na.rm = T))  # match forward primer, reverse primer, and yeast strain BCs
## [1] "matching counts (match f and r primer barcode, and moby bc): 104365740"
paste0("analysis counts: ", sum(c(poc.ce, poc.cc, faevo.ce, faevo.cc), na.rm = T))  # expected and cc counts in the analysis data
## [1] "analysis counts: 84776627"
paste0("expected analysis counts: ", sum(c(poc.ce, faevo.ce), na.rm = T))  # expected counts in the analysis data
## [1] "expected analysis counts: 81780343"
paste0("median counts per barcode: ", median(c(poc.ce, faevo.ce), na.rm = T))  # median number of counts per barcode in the expected counts dataset.
## [1] "median counts per barcode: 2961"

Analysis: Assess POC library barcode cross-contamination rate.

# calculations
# ----------------------------------------------------------------
cc = (((colSums(poc.cc, na.rm = T))/(colSums(poc.ce, na.rm = T) + colSums(poc.cc, 
    na.rm = T)))/colSums(!is.na(poc.cc))) * 100  # cc rate per exp
cc[is.nan(cc)] <- NA  # NA nans
myt <- as.data.frame(cc)  # convert to dataframe
myt$tc <- (colSums(poc.ce, na.rm = T) + colSums(poc.cc, na.rm = T))  # total counts to use as weight in weighted mean
cc001 <- weighted.mean(myt$cc, myt$tc, na.rm = T)  # save weighted mean cc (by expected reads) for grand cc rate. 
r001 <- mean(myt$tc)  # save mean for below calculation of grand cc rate. 
# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- myt$cc
myr <- myt$tc
myc <- setpalette[1]

mybinwidth <- 0.01

mytitle <- NULL
myylab <- "Count"
myxlab <- "Mean Barcode Contamination Rate (Percent Total Counts)"

myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35

mytextsize <- 8
myheight <- 3.5
mywidth <- 7

p <- ggplot(myt, aes(x = myx)) + geom_histogram(alpha = myfillapha, color = myc, 
    fill = myc, binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx, 
    myr, na.rm = T), linetype = "dashed", color = myhighlight) + geom_rug(color = myc, 
    alpha = myptalpha) + xlim(NA, 1) + theme_classic() + ggtitle(mytitle) + 
    ylab(myylab) + xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) + 
    geom_hline(yintercept = -0.01, colour = "white", size = 1.01)
p
## Warning: Removed 20 rows containing non-finite values (stat_bin).
## Warning: Removed 1 rows containing missing values (geom_bar).

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S1_A.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 20 rows containing non-finite values (stat_bin).

## Warning: Removed 1 rows containing missing values (geom_bar).
dev.off()  # 2columns wide, half as tall
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(cc, myt, p)  # remove the temporary variable

Analysis: Assess 250-generation experiment consensus library barcode cross-contamination rate.

# Calculations
# ----------------------------------------------------------------
cc = ((colSums(faevo.cc, na.rm = T))/(colSums(faevo.ce, na.rm = T) + colSums(faevo.cc, 
    na.rm = T)))/colSums(!is.na(faevo.cc)) * 100  # cc rate per exp
cc[is.nan(cc)] <- NA  # NA nans
myt <- as.data.frame(cc)  # convert to dataframe
myt$tc <- (colSums(faevo.ce, na.rm = T) + colSums(faevo.cc, na.rm = T))  # total counts to use as weight in weighted mean
cc002 <- weighted.mean(myt$cc, myt$tc, na.rm = T)  # save weighted mean cc (by expected reads) for grand cc rate. 
r002 <- mean(myt$tc)  # save mean for below calculation of grand cc rate. 

# overall cc rate
paste0("mean cross contam as % total counts = ", weighted.mean(c(cc001, cc002), 
    c(r001, r002)), " %")  # weighted mean contam weighted.mean(c(contam rate poc, contam rate faevo), c(total counts poc, total counts faevo))
## [1] "mean cross contam as % total counts = 0.0427158666120247 %"
# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- myt$cc
myr <- myt$tc
myc <- setpalette[1]

mybinwidth <- 0.01

mytitle <- NULL
myylab <- "Count"
myxlab <- "Mean Barcode Contamination Rate (Percent Total Counts)"

myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35

mytextsize <- 8
myheight <- 3.5
mywidth <- 7

p <- ggplot(myt, aes(x = myx)) + geom_histogram(alpha = myfillapha, color = myc, 
    fill = myc, binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx, 
    myr, na.rm = T), linetype = "dashed", color = myhighlight) + geom_rug(color = myc, 
    alpha = myptalpha) + xlim(NA, 1) + theme_classic() + ggtitle(mytitle) + 
    ylab(myylab) + xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) + 
    geom_hline(yintercept = -0.01, colour = "white", size = 1.01)
p
## Warning: Removed 1 rows containing missing values (geom_bar).

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S1_B.pdf", width = mywidth, height = myheight)
p
## Warning: Removed 1 rows containing missing values (geom_bar).
dev.off()  # 2columns wide, half as tall
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(cc, myt, p, cc001, cc002, r001, r002)  # remove the temporary variable

Analysis: Assess change in barcode cross-contamination rate for samples from the 250-generation experiment that were resequenced.

# Calculations & models
# -------------------------------------------------------
faevo.1.ce <- faevo.1.ce[, (colnames(faevo.1.ce) %in% colnames(faevo.2.ce))]  # subset the run 1 data to only include those columns in the rerun dataset..expected counts..
faevo.1.cc <- faevo.1.cc[, (colnames(faevo.1.cc) %in% colnames(faevo.2.cc))]  # cross contamination counts
faevo.1.ce <- faevo.1.ce[, colnames(faevo.1.ce)]
faevo.1.cc <- faevo.1.cc[, colnames(faevo.1.ce)]  # reorder the run 1 data by column name
faevo.2.ce <- faevo.2.ce[, colnames(faevo.1.ce)]
faevo.2.cc <- faevo.2.cc[, colnames(faevo.1.ce)]  # reorder the run 2 data by run 1 column name (for paired t test below)
myt1 <- ((colSums(faevo.1.cc, na.rm = T))/(colSums(faevo.1.ce, na.rm = T) + 
    colSums(faevo.1.cc, na.rm = T)))/colSums(!is.na(faevo.1.cc))  # cross contam run 1 
myt2 <- ((colSums(faevo.2.cc, na.rm = T))/(colSums(faevo.2.ce, na.rm = T) + 
    colSums(faevo.2.cc, na.rm = T)))/colSums(!is.na(faevo.2.cc))  # cross contam run 2
myt1r <- (colSums(faevo.1.ce, na.rm = T) + colSums(faevo.1.cc, na.rm = T))  # reads run 1
myt2r <- (colSums(faevo.2.ce, na.rm = T) + colSums(faevo.2.cc, na.rm = T))  # reads run 2

# unweighted t-test results match weighted t-test qualitatively. see
# weighted t.test imemediately below.  var.test(myt1, myt2) # variance is
# not equal mymod <- t.test(myt1, myt2, alternative = 'greater', paired =
# TRUE, var.equal = FALSE); mymod # t.test, paired, 1 sided (greater),
# unequal variance

myt <- myt1 - myt2  # calculatation for weighted t.test below (contam run 1 - contam run 2)
mytr <- 2/((1/myt1r) + (1/myt2r))  # reads for calculation below (harmonic mean of reads for run 1 and run 2)
mymod2 <- wtd.t.test(x = myt, y = 0, weight = mytr, alternative = "greater", 
    mean1 = T)  # t.test, paired, 1 sided (greater), weighted
# -----------------------------------------------------------------------------

# create temporary plotting dataframe
# -----------------------------------------
myt1 <- as.data.frame(cbind(myt1, myt1r))
colnames(myt1) <- c("cc", "r")
myt1$seq.run <- 1  # run 1 df
myt2 <- as.data.frame(cbind(myt2, myt2r))
colnames(myt2) <- c("cc", "r")
myt2$seq.run <- 2  # run 2 df
myt <- rbind(myt1, myt2)
myt$seq.run <- as.factor(myt$seq.run)  # combine
# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- myt$seq.run
myy <- myt$cc
myr <- myt$r

myx1 <- myt1$cc
myr1 <- myt1$r

myx2 <- myt2$cc
myr2 <- myt2$r

mytitle <- NULL
myylab <- "Mean Barcode Contamination Rate (% Total Counts)"
myxlab <- "Sequenced Libary Run Number"

myoutline <- "grey50"
myfill <- "grey90"
myhighlight <- "grey50"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette

mytextsize <- 8
myheight <- 7
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(myt, aes(x = myx, y = myy, weight = myr, color = myx)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_segment(x = 0.75, xend = 1.25, y = weighted.mean(myx1, 
    myr1), yend = weighted.mean(myx1, myr1), color = myhighlight) + geom_segment(x = 1.75, 
    xend = 2.25, y = weighted.mean(myx2, myr2), yend = weighted.mean(myx2, myr2), 
    color = myhighlight) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, 
    size = myptsize) + geom_rug(sides = "l", alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + theme_classic() + theme(legend.position = "none") + 
    ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + scale_x_discrete(breaks = c("1", 
    "2"), labels = c("Library 1", "Library 2")) + theme(axis.title.y = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0))
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S2.pdf", width = mywidth, height = myheight)
p
dev.off()  # 7by7
## quartz_off_screen 
##                 2

# -----------------------------------------------------------------------------
rm(faevo.1.ce, faevo.1.cc, faevo.2.ce, faevo.2.cc, myt1, myt1r, myt2, myt2r, 
    myt, mytr, p)  # clean up this section. 
rm(faevo.cc, faevo.ce, poc.cc, poc.ce)  # clean up for next major set of analyses (remove the counts datasets to keep our global environment clean)

Analysis: Assess whether detected fitness change for barcodes in the 250-generation experiment is affected by barcode-cross contamination rate

my <- myfa[order(myfa$bctID), ]
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lmer(scale(deltafit) ~ scale(ccdeltafit) + (1 | treatment/bctID), weights = rdeltafit, 
    data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
summary(mymod)  # problems with convergence here...solved with an optimizer and some extra iterations.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(deltafit) ~ scale(ccdeltafit) + (1 | treatment/bctID)
##    Data: my
## Weights: rdeltafit
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
## 
## REML criterion at convergence: 797.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1167 -0.4622  0.0163  0.6366  3.5030 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  bctID:treatment (Intercept)   0.1174  0.3426 
##  treatment       (Intercept)   0.1865  0.4319 
##  Residual                    368.2615 19.1901 
## Number of obs: 608, groups:  bctID:treatment, 152; treatment, 6
## 
## Fixed effects:
##                    Estimate Std. Error        df t value
## (Intercept)        -0.11495    0.17940   4.94812  -0.641
## scale(ccdeltafit)  -0.26885    0.01858 579.59718 -14.473
##                              Pr(>|t|)    
## (Intercept)                      0.55    
## scale(ccdeltafit) <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## scl(ccdltf) 0.007
# mymod <- lm(deltafit ~ ccdeltafit, weights = rdeltafit, data =
# my);r.squaredGLMM(mymod) # get the r2 for the lm version of the model if
# desired.
hist(residuals(mymod), breaks = 100)  # residuals have outlying points -- but biew the weighted version of the residual histogram as well....

hist(expand_residuals(mymod, my$rdeltafit), breaks = 100)  # this more accurately represents the distribution of residuals in the model output.

# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    show.r2 = F, file = "000_Table_S6.html")

	scale(deltafit)
Predictors	Estimates	CI	Statistic	p
(Intercept)	-0.11495	-0.46656 – 0.23667	-0.64073	0.55018888
scale(ccdeltafit)	-0.26885	-0.30526 – -0.23244	-14.47319	<0.001
Random Effects
σ2	368.26
τ00 bctID:treatment	0.12
τ00 treatment	0.19
ICC	0.00
N bctID	152
N treatment	6
Observations	608

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
mymod <- lmer(deltafit ~ ccdeltafit + (1 | treatment/bctID), weights = rdeltafit, 
    data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))  # unmanipulated data for plotting
## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling
myx <- my$ccdeltafit
myy <- my$deltafit
myr <- my$rdeltafit
myc <- as.character(my$treatment)
myc[myc == "CM+Ethanoldiploid0.001"] <- "CM+EtOH \n Diploid \n 0.001 \n"
myc[myc == "CM+Saltdiploid0.001"] <- "CM+NaCl \n Diploid \n 0.001 \n"
myc[myc == "CMdiploid0.001"] <- "CM \n Diploid \n 0.001 \n"
myc[myc == "CMdiploid0.004"] <- "CM \n Diploid \n 0.004 \n"
myc[myc == "CMdiploid0.00025"] <- "CM \n Diploid \n 0.00025 \n"
myc[myc == "CMhaploid0.001"] <- "CM \n Haploid \n 0.001 \n"
myc <- factor(myc, levels = c("CM \n Diploid \n 0.001 \n", "CM+EtOH \n Diploid \n 0.001 \n", 
    "CM+NaCl \n Diploid \n 0.001 \n", "CM \n Diploid \n 0.00025 \n", "CM \n Diploid \n 0.004 \n", 
    "CM \n Haploid \n 0.001 \n"))

myg <- paste0(my$treatment, my$bctID)

mytitle <- NULL
myylab <- "Change in Fitness"
myxlab <- "Mean Barcode Contamination Rate (% Total Counts)"

myhighlight <- "black"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette

mytextsize <- 8
myheight <- 7
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, color = myc)) + geom_line(aes(y = predict(mymod), 
    group = myg, color = myc), size = 1, alpha = myptalpha) + geom_point(aes(color = myc), 
    alpha = myptalpha, size = myptsize) + geom_rug(alpha = myptalpha) + scale_color_manual(values = setpalette) + 
    theme_classic() + labs(color = "Treatment") + ggtitle(mytitle) + ylab(myylab) + 
    xlab(myxlab) + theme(axis.title.y = element_text(size = mytextsize, face = "bold", 
    margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) + 
    theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S3.pdf", width = mywidth, height = myheight)
p
dev.off()  # 7by7
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p, weights)

Analysis: Assess fitness change in 250-generations of experimental evolution for barcodes in the 250-generation experiment.

my <- myfa[order(myfa$bctID), ]
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lm(deltafit ~ bctID + 0, weights = rdeltafit, data = my)
summary(mymod)
## 
## Call:
## lm(formula = deltafit ~ bctID + 0, data = my, weights = rdeltafit)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.9024 -0.9879  0.0734  0.9561  4.8229 
## 
## Coefficients:
##                                     Estimate  Std. Error t value
## bctIDCM+Ethanoldiploid0.001d1A2   0.00224666  0.01223756   0.184
## bctIDCM+Ethanoldiploid0.001d1A4   0.01765500  0.01365937   1.293
## bctIDCM+Ethanoldiploid0.001d1A5  -0.00183138  0.01601185  -0.114
## bctIDCM+Ethanoldiploid0.001d1A6   0.00847721  0.01160185   0.731
## bctIDCM+Ethanoldiploid0.001d1A7   0.00837054  0.01577339   0.531
## bctIDCM+Ethanoldiploid0.001d1D10  0.03354979  0.01054805   3.181
## bctIDCM+Ethanoldiploid0.001d1D11  0.02980938  0.01093362   2.726
## bctIDCM+Ethanoldiploid0.001d1D12  0.02399057  0.01018367   2.356
## bctIDCM+Ethanoldiploid0.001d1D9   0.02343344  0.02027118   1.156
## bctIDCM+Ethanoldiploid0.001d1E10  0.04807661  0.01120179   4.292
## bctIDCM+Ethanoldiploid0.001d1E11  0.03162314  0.01125755   2.809
## bctIDCM+Ethanoldiploid0.001d1E12 -0.03049231  0.01211023  -2.518
## bctIDCM+Ethanoldiploid0.001d1E9  -0.00999714  0.01070592  -0.934
## bctIDCM+Ethanoldiploid0.001d1F1   0.00938513  0.00879011   1.068
## bctIDCM+Ethanoldiploid0.001d1F2   0.01938389  0.01121524   1.728
## bctIDCM+Ethanoldiploid0.001d1G1   0.01657615  0.01500470   1.105
## bctIDCM+Ethanoldiploid0.001d1G2  -0.01127658  0.01818702  -0.620
## bctIDCM+Ethanoldiploid0.001d1H2  -0.01386016  0.01837215  -0.754
## bctIDCM+Ethanoldiploid0.001d1H4   0.00296372  0.02376768   0.125
## bctIDCM+Ethanoldiploid0.001d1H5   0.00347374  0.01381229   0.251
## bctIDCM+Ethanoldiploid0.001d1H6   0.00333766  0.01279047   0.261
## bctIDCM+Ethanoldiploid0.001d1H7  -0.00551949  0.02558074  -0.216
## bctIDCM+Saltdiploid0.001d1A2      0.04403600  0.01205900   3.652
## bctIDCM+Saltdiploid0.001d1A3      0.01686781  0.01180626   1.429
## bctIDCM+Saltdiploid0.001d1A4      0.02820342  0.01240104   2.274
## bctIDCM+Saltdiploid0.001d1A5      0.39852910  0.19667177   2.026
## bctIDCM+Saltdiploid0.001d1A7      0.02578701  0.01322947   1.949
## bctIDCM+Saltdiploid0.001d1F3      0.14079891  0.02033980   6.922
## bctIDCM+Saltdiploid0.001d1F4      0.18806044  0.02296431   8.189
## bctIDCM+Saltdiploid0.001d1F5      0.19913858  0.15799739   1.260
## bctIDCM+Saltdiploid0.001d1F6      0.23451078  0.01480773  15.837
## bctIDCM+Saltdiploid0.001d1F7      0.15255139  0.02121079   7.192
## bctIDCM+Saltdiploid0.001d1F8      0.12888704  0.01760927   7.319
## bctIDCM+Saltdiploid0.001d1G3      0.10101471  0.03918185   2.578
## bctIDCM+Saltdiploid0.001d1G4      0.13402619  0.04517263   2.967
## bctIDCM+Saltdiploid0.001d1G5      0.09585476  0.01936410   4.950
## bctIDCM+Saltdiploid0.001d1G6      0.08971402  0.03285475   2.731
## bctIDCM+Saltdiploid0.001d1G7      0.08740043  0.02106254   4.150
## bctIDCM+Saltdiploid0.001d1G8      0.12572385  0.03290001   3.821
## bctIDCM+Saltdiploid0.001d1H2      0.00495816  0.01539561   0.322
## bctIDCM+Saltdiploid0.001d1H3      0.02115633  0.01617586   1.308
## bctIDCM+Saltdiploid0.001d1H4      0.02137396  0.01932366   1.106
## bctIDCM+Saltdiploid0.001d1H5      0.00354951  0.01200179   0.296
## bctIDCM+Saltdiploid0.001d1H7      0.09417670  0.02624046   3.589
## bctIDCMdiploid0.00025d1A2         0.01559982  0.01599335   0.975
## bctIDCMdiploid0.00025d1A3        -0.00878718  0.01495776  -0.587
## bctIDCMdiploid0.00025d1A4         0.00174608  0.01441827   0.121
## bctIDCMdiploid0.00025d1A5         0.02677688  0.01511949   1.771
## bctIDCMdiploid0.00025d1A6         0.01198957  0.01303287   0.920
## bctIDCMdiploid0.00025d1A7        -0.00440271  0.01621492  -0.272
## bctIDCMdiploid0.00025d1B10       -0.00809898  0.01653299  -0.490
## bctIDCMdiploid0.00025d1B11        0.01723177  0.01242047   1.387
## bctIDCMdiploid0.00025d1B12        0.02983743  0.01830367   1.630
## bctIDCMdiploid0.00025d1B9         0.02155039  0.01336058   1.613
## bctIDCMdiploid0.00025d1C10       -0.01674094  0.01859239  -0.900
## bctIDCMdiploid0.00025d1C11        0.03216629  0.01622343   1.983
## bctIDCMdiploid0.00025d1C12        0.01496869  0.01326625   1.128
## bctIDCMdiploid0.00025d1C9        -0.00518047  0.01188006  -0.436
## bctIDCMdiploid0.00025d1D2         0.00370897  0.05550949   0.067
## bctIDCMdiploid0.00025d1E2         0.01289359  0.01052145   1.225
## bctIDCMdiploid0.00025d1H2         0.00760594  0.01762321   0.432
## bctIDCMdiploid0.00025d1H3        -0.00458899  0.01745253  -0.263
## bctIDCMdiploid0.00025d1H4         0.01843903  0.02426240   0.760
## bctIDCMdiploid0.00025d1H5         0.02541993  0.01377634   1.845
## bctIDCMdiploid0.00025d1H6         0.00414724  0.01352400   0.307
## bctIDCMdiploid0.00025d1H7         0.00190764  0.02782539   0.069
## bctIDCMdiploid0.001d1A11          0.03612225  0.01815435   1.990
## bctIDCMdiploid0.001d1A2           0.00809289  0.01417709   0.571
## bctIDCMdiploid0.001d1A2B          0.00927271  0.01214216   0.764
## bctIDCMdiploid0.001d1A3           0.00032174  0.01485237   0.022
## bctIDCMdiploid0.001d1A3B          0.00398055  0.01385717   0.287
## bctIDCMdiploid0.001d1A4           0.01780845  0.01310805   1.359
## bctIDCMdiploid0.001d1A4B          0.00153330  0.01338588   0.115
## bctIDCMdiploid0.001d1A5           0.02937872  0.01569201   1.872
## bctIDCMdiploid0.001d1A5B         -0.00512256  0.01465871  -0.349
## bctIDCMdiploid0.001d1A6           0.02327336  0.01302785   1.786
## bctIDCMdiploid0.001d1A6B          0.00955064  0.01339667   0.713
## bctIDCMdiploid0.001d1A7           0.04290593  0.01344789   3.191
## bctIDCMdiploid0.001d1A8           0.01512575  0.01169949   1.293
## bctIDCMdiploid0.001d1A9           0.02268337  0.01128163   2.011
## bctIDCMdiploid0.001d1B1           0.05232713  0.01916650   2.730
## bctIDCMdiploid0.001d1B2           0.04063255  0.01363624   2.980
## bctIDCMdiploid0.001d1C1           0.05496896  0.01813523   3.031
## bctIDCMdiploid0.001d1C2           0.03217244  0.01299148   2.476
## bctIDCMdiploid0.001d1D4           0.08813997  0.05292837   1.665
## bctIDCMdiploid0.001d1D5          -0.01879019  0.01310839  -1.433
## bctIDCMdiploid0.001d1D6           0.04266468  0.04243456   1.005
## bctIDCMdiploid0.001d1D7           0.01221091  0.01202676   1.015
## bctIDCMdiploid0.001d1D8          -0.01270792  0.01081313  -1.175
## bctIDCMdiploid0.001d1E4           0.00301728  0.00782870   0.385
## bctIDCMdiploid0.001d1E5           0.03661495  0.02836404   1.291
## bctIDCMdiploid0.001d1E6          -0.02094510  0.00919326  -2.278
## bctIDCMdiploid0.001d1E7           0.02727065  0.01036588   2.631
## bctIDCMdiploid0.001d1E8          -0.00435115  0.01093102  -0.398
## bctIDCMdiploid0.001d1H11          0.03103102  0.01491899   2.080
## bctIDCMdiploid0.001d1H2           0.01207171  0.01797336   0.672
## bctIDCMdiploid0.001d1H2B         -0.01771525  0.01722331  -1.029
## bctIDCMdiploid0.001d1H3          -0.00572467  0.01932551  -0.296
## bctIDCMdiploid0.001d1H3B          0.00004411  0.01824231   0.002
## bctIDCMdiploid0.001d1H4           0.00887392  0.02490203   0.356
## bctIDCMdiploid0.001d1H4B         -0.00268735  0.02541084  -0.106
## bctIDCMdiploid0.001d1H5           0.01359003  0.01528546   0.889
## bctIDCMdiploid0.001d1H5B         -0.00274906  0.01359686  -0.202
## bctIDCMdiploid0.001d1H6           0.03458717  0.01367055   2.530
## bctIDCMdiploid0.001d1H6B          0.02557412  0.01329414   1.924
## bctIDCMdiploid0.001d1H7           0.00942620  0.02647801   0.356
## bctIDCMdiploid0.001d1H8           0.01372130  0.02208970   0.621
## bctIDCMdiploid0.001d1H9           0.01922550  0.02317816   0.829
## bctIDCMdiploid0.004d1A2           0.00073845  0.01697648   0.043
## bctIDCMdiploid0.004d1A3           0.01219203  0.01903325   0.641
## bctIDCMdiploid0.004d1A4           0.05056340  0.01731766   2.920
## bctIDCMdiploid0.004d1A5           0.05149485  0.01798080   2.864
## bctIDCMdiploid0.004d1A6           0.01931859  0.01722893   1.121
## bctIDCMdiploid0.004d1A7           0.02865780  0.02004130   1.430
## bctIDCMdiploid0.004d1B3           0.02462782  0.01445493   1.704
## bctIDCMdiploid0.004d1B4           0.03125331  0.01659619   1.883
## bctIDCMdiploid0.004d1B5           0.04298949  0.01367318   3.144
## bctIDCMdiploid0.004d1B6           0.05327518  0.01738231   3.065
## bctIDCMdiploid0.004d1B7           0.07749931  0.01526527   5.077
## bctIDCMdiploid0.004d1C3           0.04273794  0.01469507   2.908
## bctIDCMdiploid0.004d1C4           0.02312320  0.01388803   1.665
## bctIDCMdiploid0.004d1C5           0.04334705  0.01293091   3.352
## bctIDCMdiploid0.004d1C6           0.06362223  0.01397915   4.551
## bctIDCMdiploid0.004d1C7           0.03835768  0.02896918   1.324
## bctIDCMdiploid0.004d1H2           0.02225147  0.02308198   0.964
## bctIDCMdiploid0.004d1H3           0.02106572  0.02173722   0.969
## bctIDCMdiploid0.004d1H4           0.00548235  0.03044485   0.180
## bctIDCMdiploid0.004d1H5           0.06026645  0.01911119   3.153
## bctIDCMdiploid0.004d1H6           0.04022208  0.01766259   2.277
## bctIDCMdiploid0.004d1H7           0.01164662  0.03155608   0.369
## bctIDCMhaploid0.001d1A2           0.17748276  0.18501465   0.959
## bctIDCMhaploid0.001d1A3           0.03475922  0.01642162   2.117
## bctIDCMhaploid0.001d1A3B          0.02653209  0.01137336   2.333
## bctIDCMhaploid0.001d1A4           0.07888998  0.01157424   6.816
## bctIDCMhaploid0.001d1A4B          0.00800659  0.01199416   0.668
## bctIDCMhaploid0.001d1A5           0.07567154  0.01371413   5.518
## bctIDCMhaploid0.001d1A5B          0.05376713  0.01259551   4.269
## bctIDCMhaploid0.001d1A6           0.02913551  0.01052784   2.767
## bctIDCMhaploid0.001d1A6B          0.02748147  0.00908708   3.024
## bctIDCMhaploid0.001d1A7           0.28424691  0.22971593   1.237
## bctIDCMhaploid0.001d1A7B          0.37549211  0.27832375   1.349
## bctIDCMhaploid0.001d1H2           0.05734718  0.01492184   3.843
## bctIDCMhaploid0.001d1H3           0.00588798  0.01683299   0.350
## bctIDCMhaploid0.001d1H3B          0.02271945  0.01846922   1.230
## bctIDCMhaploid0.001d1H4           0.06491278  0.03830520   1.695
## bctIDCMhaploid0.001d1H4B          0.01454690  0.01990335   0.731
## bctIDCMhaploid0.001d1H5           0.24925362  0.21028584   1.185
## bctIDCMhaploid0.001d1H5B          0.07768822  0.01659155   4.682
## bctIDCMhaploid0.001d1H6           0.10010351  0.02954796   3.388
## bctIDCMhaploid0.001d1H6B          0.11256605  0.02312772   4.867
## bctIDCMhaploid0.001d1H7           0.08865934  0.02084663   4.253
## bctIDCMhaploid0.001d1H7B          0.07939256  0.01932983   4.107
##                                              Pr(>|t|)    
## bctIDCM+Ethanoldiploid0.001d1A2              0.854419    
## bctIDCM+Ethanoldiploid0.001d1A4              0.196832    
## bctIDCM+Ethanoldiploid0.001d1A5              0.908990    
## bctIDCM+Ethanoldiploid0.001d1A6              0.465352    
## bctIDCM+Ethanoldiploid0.001d1A7              0.595903    
## bctIDCM+Ethanoldiploid0.001d1D10             0.001570 ** 
## bctIDCM+Ethanoldiploid0.001d1D11             0.006650 ** 
## bctIDCM+Ethanoldiploid0.001d1D12             0.018906 *  
## bctIDCM+Ethanoldiploid0.001d1D9              0.248287    
## bctIDCM+Ethanoldiploid0.001d1E10  0.00002165288566661 ***
## bctIDCM+Ethanoldiploid0.001d1E11             0.005182 ** 
## bctIDCM+Ethanoldiploid0.001d1E12             0.012147 *  
## bctIDCM+Ethanoldiploid0.001d1E9              0.350903    
## bctIDCM+Ethanoldiploid0.001d1F1              0.286224    
## bctIDCM+Ethanoldiploid0.001d1F2              0.084602 .  
## bctIDCM+Ethanoldiploid0.001d1G1              0.269859    
## bctIDCM+Ethanoldiploid0.001d1G2              0.535545    
## bctIDCM+Ethanoldiploid0.001d1H2              0.450992    
## bctIDCM+Ethanoldiploid0.001d1H4              0.900820    
## bctIDCM+Ethanoldiploid0.001d1H5              0.801543    
## bctIDCM+Ethanoldiploid0.001d1H6              0.794250    
## bctIDCM+Ethanoldiploid0.001d1H7              0.829265    
## bctIDCM+Saltdiploid0.001d1A2                 0.000291 ***
## bctIDCM+Saltdiploid0.001d1A3                 0.153770    
## bctIDCM+Saltdiploid0.001d1A4                 0.023413 *  
## bctIDCM+Saltdiploid0.001d1A5                 0.043309 *  
## bctIDCM+Saltdiploid0.001d1A7                 0.051883 .  
## bctIDCM+Saltdiploid0.001d1F3      0.00000000001514729 ***
## bctIDCM+Saltdiploid0.001d1F4      0.00000000000000266 ***
## bctIDCM+Saltdiploid0.001d1F5                 0.208173    
## bctIDCM+Saltdiploid0.001d1F6     < 0.0000000000000002 ***
## bctIDCM+Saltdiploid0.001d1F7      0.00000000000263293 ***
## bctIDCM+Saltdiploid0.001d1F8      0.00000000000113446 ***
## bctIDCM+Saltdiploid0.001d1G3                 0.010247 *  
## bctIDCM+Saltdiploid0.001d1G4                 0.003166 ** 
## bctIDCM+Saltdiploid0.001d1G5      0.00000104573824082 ***
## bctIDCM+Saltdiploid0.001d1G6                 0.006566 ** 
## bctIDCM+Saltdiploid0.001d1G7      0.00003974952106726 ***
## bctIDCM+Saltdiploid0.001d1G8                 0.000151 ***
## bctIDCM+Saltdiploid0.001d1H2                 0.747562    
## bctIDCM+Saltdiploid0.001d1H3                 0.191568    
## bctIDCM+Saltdiploid0.001d1H4                 0.269265    
## bctIDCM+Saltdiploid0.001d1H5                 0.767557    
## bctIDCM+Saltdiploid0.001d1H7                 0.000368 ***
## bctIDCMdiploid0.00025d1A2                    0.329882    
## bctIDCMdiploid0.00025d1A3                    0.557181    
## bctIDCMdiploid0.00025d1A4                    0.903664    
## bctIDCMdiploid0.00025d1A5                    0.077226 .  
## bctIDCMdiploid0.00025d1A6                    0.358086    
## bctIDCMdiploid0.00025d1A7                    0.786112    
## bctIDCMdiploid0.00025d1B10                   0.624463    
## bctIDCMdiploid0.00025d1B11                   0.166007    
## bctIDCMdiploid0.00025d1B12                   0.103764    
## bctIDCMdiploid0.00025d1B9                    0.107440    
## bctIDCMdiploid0.00025d1C10                   0.368373    
## bctIDCMdiploid0.00025d1C11                   0.048000 *  
## bctIDCMdiploid0.00025d1C12                   0.259774    
## bctIDCMdiploid0.00025d1C9                    0.662997    
## bctIDCMdiploid0.00025d1D2                    0.946757    
## bctIDCMdiploid0.00025d1E2                    0.221035    
## bctIDCMdiploid0.00025d1H2                    0.666246    
## bctIDCMdiploid0.00025d1H3                    0.792714    
## bctIDCMdiploid0.00025d1H4                    0.447657    
## bctIDCMdiploid0.00025d1H5                    0.065659 .  
## bctIDCMdiploid0.00025d1H6                    0.759244    
## bctIDCMdiploid0.00025d1H7                    0.945372    
## bctIDCMdiploid0.001d1A11                     0.047218 *  
## bctIDCMdiploid0.001d1A2                      0.568387    
## bctIDCMdiploid0.001d1A2B                     0.445453    
## bctIDCMdiploid0.001d1A3                      0.982727    
## bctIDCMdiploid0.001d1A3B                     0.774047    
## bctIDCMdiploid0.001d1A4                      0.174949    
## bctIDCMdiploid0.001d1A4B                     0.908856    
## bctIDCMdiploid0.001d1A5                      0.061817 .  
## bctIDCMdiploid0.001d1A5B                     0.726909    
## bctIDCMdiploid0.001d1A6                      0.074694 .  
## bctIDCMdiploid0.001d1A6B                     0.476265    
## bctIDCMdiploid0.001d1A7                      0.001518 ** 
## bctIDCMdiploid0.001d1A8                      0.196716    
## bctIDCMdiploid0.001d1A9                      0.044951 *  
## bctIDCMdiploid0.001d1B1                      0.006576 ** 
## bctIDCMdiploid0.001d1B2                      0.003039 ** 
## bctIDCMdiploid0.001d1C1                      0.002576 ** 
## bctIDCMdiploid0.001d1C2                      0.013632 *  
## bctIDCMdiploid0.001d1D4                      0.096546 .  
## bctIDCMdiploid0.001d1D5                      0.152415    
## bctIDCMdiploid0.001d1D6                      0.315226    
## bctIDCMdiploid0.001d1D7                      0.310496    
## bctIDCMdiploid0.001d1D8                      0.240516    
## bctIDCMdiploid0.001d1E4                      0.700111    
## bctIDCMdiploid0.001d1E5                      0.197395    
## bctIDCMdiploid0.001d1E6                      0.023170 *  
## bctIDCMdiploid0.001d1E7                      0.008807 ** 
## bctIDCMdiploid0.001d1E8                      0.690775    
## bctIDCMdiploid0.001d1H11                     0.038086 *  
## bctIDCMdiploid0.001d1H2                      0.502150    
## bctIDCMdiploid0.001d1H2B                     0.304231    
## bctIDCMdiploid0.001d1H3                      0.767194    
## bctIDCMdiploid0.001d1H3B                     0.998072    
## bctIDCMdiploid0.001d1H4                      0.721741    
## bctIDCMdiploid0.001d1H4B                     0.915822    
## bctIDCMdiploid0.001d1H5                      0.374428    
## bctIDCMdiploid0.001d1H5B                     0.839864    
## bctIDCMdiploid0.001d1H6                      0.011740 *  
## bctIDCMdiploid0.001d1H6B                     0.055013 .  
## bctIDCMdiploid0.001d1H7                      0.722004    
## bctIDCMdiploid0.001d1H8                      0.534803    
## bctIDCMdiploid0.001d1H9                      0.407275    
## bctIDCMdiploid0.004d1A2                      0.965323    
## bctIDCMdiploid0.004d1A3                      0.522127    
## bctIDCMdiploid0.004d1A4                      0.003677 ** 
## bctIDCMdiploid0.004d1A5                      0.004378 ** 
## bctIDCMdiploid0.004d1A6                      0.262755    
## bctIDCMdiploid0.004d1A7                      0.153420    
## bctIDCMdiploid0.004d1B3                      0.089106 .  
## bctIDCMdiploid0.004d1B4                      0.060315 .  
## bctIDCMdiploid0.004d1B5                      0.001775 ** 
## bctIDCMdiploid0.004d1B6                      0.002307 ** 
## bctIDCMdiploid0.004d1B7           0.00000056011850773 ***
## bctIDCMdiploid0.004d1C3                      0.003811 ** 
## bctIDCMdiploid0.004d1C4                      0.096605 .  
## bctIDCMdiploid0.004d1C5                      0.000868 ***
## bctIDCMdiploid0.004d1C6           0.00000684943418923 ***
## bctIDCMdiploid0.004d1C7                      0.186138    
## bctIDCMdiploid0.004d1H2                      0.335547    
## bctIDCMdiploid0.004d1H3                      0.333005    
## bctIDCMdiploid0.004d1H4                      0.857174    
## bctIDCMdiploid0.004d1H5                      0.001720 ** 
## bctIDCMdiploid0.004d1H6                      0.023234 *  
## bctIDCMdiploid0.004d1H7                      0.712242    
## bctIDCMhaploid0.001d1A2                      0.337921    
## bctIDCMhaploid0.001d1A3                      0.034829 *  
## bctIDCMhaploid0.001d1A3B                     0.020091 *  
## bctIDCMhaploid0.001d1A4           0.00000000002976149 ***
## bctIDCMhaploid0.001d1A4B                     0.504764    
## bctIDCMhaploid0.001d1A5           0.00000005766073013 ***
## bctIDCMhaploid0.001d1A5B          0.00002392620167394 ***
## bctIDCMhaploid0.001d1A6                      0.005879 ** 
## bctIDCMhaploid0.001d1A6B                     0.002633 ** 
## bctIDCMhaploid0.001d1A7                      0.216581    
## bctIDCMhaploid0.001d1A7B                     0.177968    
## bctIDCMhaploid0.001d1H2                      0.000139 ***
## bctIDCMhaploid0.001d1H3                      0.726659    
## bctIDCMhaploid0.001d1H3B                     0.219285    
## bctIDCMhaploid0.001d1H4                      0.090830 .  
## bctIDCMhaploid0.001d1H4B                     0.465230    
## bctIDCMhaploid0.001d1H5                      0.236513    
## bctIDCMhaploid0.001d1H5B          0.00000374585588773 ***
## bctIDCMhaploid0.001d1H6                      0.000766 ***
## bctIDCMhaploid0.001d1H6B          0.00000156280464386 ***
## bctIDCMhaploid0.001d1H7           0.00002561158235811 ***
## bctIDCMhaploid0.001d1H7B          0.00004746151175424 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.922 on 456 degrees of freedom
## Multiple R-squared:  0.7186, Adjusted R-squared:  0.6248 
## F-statistic:  7.66 on 152 and 456 DF,  p-value: < 0.00000000000000022
hist(residuals(mymod), breaks = 100)  # tails do look a little heavy

hist(expand_residuals(mymod, my$rdeltafit), breaks = 100)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # model fit (multi) -- looks pretty good. 

# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S7.html")  # table does not have corrected p-values. 

	deltafit
Predictors	Estimates	CI	Statistic	p
CM+Ethanoldiploid 0.001 d 1 A 2	0.00225	-0.02174 – 0.02623	0.18359	0.85441888
CM+Ethanoldiploid 0.001 d 1 A 4	0.01765	-0.00912 – 0.04443	1.29252	0.19683180
CM+Ethanoldiploid 0.001 d 1 A 5	-0.00183	-0.03321 – 0.02955	-0.11438	0.90898955
CM+Ethanoldiploid 0.001 d 1 A 6	0.00848	-0.01426 – 0.03122	0.73068	0.46535168
CM+Ethanoldiploid 0.001 d 1 A 7	0.00837	-0.02254 – 0.03929	0.53067	0.59590266
CM+Ethanoldiploid 0.001 d 1 D 10	0.03355	0.01288 – 0.05422	3.18066	0.00156959
CM+Ethanoldiploid 0.001 d 1 D 11	0.02981	0.00838 – 0.05124	2.72640	0.00664970
CM+Ethanoldiploid 0.001 d 1 D 12	0.02399	0.00403 – 0.04395	2.35579	0.01890587
CM+Ethanoldiploid 0.001 d 1 D 9	0.02343	-0.01630 – 0.06316	1.15600	0.24828733
CM+Ethanoldiploid 0.001 d 1 E 10	0.04808	0.02612 – 0.07003	4.29187	0.00002165
CM+Ethanoldiploid 0.001 d 1 E 11	0.03162	0.00956 – 0.05369	2.80906	0.00518210
CM+Ethanoldiploid 0.001 d 1 E 12	-0.03049	-0.05423 – -0.00676	-2.51790	0.01214729
CM+Ethanoldiploid 0.001 d 1 E 9	-0.01000	-0.03098 – 0.01099	-0.93380	0.35090332
CM+Ethanoldiploid 0.001 d 1 F 1	0.00939	-0.00784 – 0.02661	1.06769	0.28622445
CM+Ethanoldiploid 0.001 d 1 F 2	0.01938	-0.00260 – 0.04137	1.72835	0.08460206
CM+Ethanoldiploid 0.001 d 1 G 1	0.01658	-0.01283 – 0.04598	1.10473	0.26985901
CM+Ethanoldiploid 0.001 d 1 G 2	-0.01128	-0.04692 – 0.02437	-0.62003	0.53554478
CM+Ethanoldiploid 0.001 d 1 H 2	-0.01386	-0.04987 – 0.02215	-0.75441	0.45099161
CM+Ethanoldiploid 0.001 d 1 H 4	0.00296	-0.04362 – 0.04955	0.12470	0.90081956
CM+Ethanoldiploid 0.001 d 1 H 5	0.00347	-0.02360 – 0.03055	0.25150	0.80154341
CM+Ethanoldiploid 0.001 d 1 H 6	0.00334	-0.02173 – 0.02841	0.26095	0.79424982
CM+Ethanoldiploid 0.001 d 1 H 7	-0.00552	-0.05566 – 0.04462	-0.21577	0.82926550
CM+Saltdiploid 0.001 d 1 A 2	0.04404	0.02040 – 0.06767	3.65171	0.00029070
CM+Saltdiploid 0.001 d 1 A 3	0.01687	-0.00627 – 0.04001	1.42872	0.15376998
CM+Saltdiploid 0.001 d 1 A 4	0.02820	0.00390 – 0.05251	2.27428	0.02341305
CM+Saltdiploid 0.001 d 1 A 5	0.39853	0.01306 – 0.78400	2.02637	0.04330880
CM+Saltdiploid 0.001 d 1 A 7	0.02579	-0.00014 – 0.05172	1.94921	0.05188310
CM+Saltdiploid 0.001 d 1 F 3	0.14080	0.10093 – 0.18066	6.92234	<0.001
CM+Saltdiploid 0.001 d 1 F 4	0.18806	0.14305 – 0.23307	8.18925	<0.001
CM+Saltdiploid 0.001 d 1 F 5	0.19914	-0.11053 – 0.50881	1.26039	0.20817273
CM+Saltdiploid 0.001 d 1 F 6	0.23451	0.20549 – 0.26353	15.83705	<0.001
CM+Saltdiploid 0.001 d 1 F 7	0.15255	0.11098 – 0.19412	7.19216	<0.001
CM+Saltdiploid 0.001 d 1 F 8	0.12889	0.09437 – 0.16340	7.31927	<0.001
CM+Saltdiploid 0.001 d 1 G 3	0.10101	0.02422 – 0.17781	2.57810	0.01024726
CM+Saltdiploid 0.001 d 1 G 4	0.13403	0.04549 – 0.22256	2.96698	0.00316552
CM+Saltdiploid 0.001 d 1 G 5	0.09585	0.05790 – 0.13381	4.95013	0.00000105
CM+Saltdiploid 0.001 d 1 G 6	0.08971	0.02532 – 0.15411	2.73063	0.00656636
CM+Saltdiploid 0.001 d 1 G 7	0.08740	0.04612 – 0.12868	4.14957	0.00003975
CM+Saltdiploid 0.001 d 1 G 8	0.12572	0.06124 – 0.19021	3.82139	0.00015111
CM+Saltdiploid 0.001 d 1 H 2	0.00496	-0.02522 – 0.03513	0.32205	0.74756239
CM+Saltdiploid 0.001 d 1 H 3	0.02116	-0.01055 – 0.05286	1.30790	0.19156762
CM+Saltdiploid 0.001 d 1 H 4	0.02137	-0.01650 – 0.05925	1.10610	0.26926505
CM+Saltdiploid 0.001 d 1 H 5	0.00355	-0.01997 – 0.02707	0.29575	0.76755679
CM+Saltdiploid 0.001 d 1 H 7	0.09418	0.04275 – 0.14561	3.58899	0.00036794
C Mdiploid 0.00025 d 1 A 2	0.01560	-0.01575 – 0.04695	0.97539	0.32988206
C Mdiploid 0.00025 d 1 A 3	-0.00879	-0.03810 – 0.02053	-0.58747	0.55718131
C Mdiploid 0.00025 d 1 A 4	0.00175	-0.02651 – 0.03001	0.12110	0.90366385
C Mdiploid 0.00025 d 1 A 5	0.02678	-0.00286 – 0.05641	1.77102	0.07722576
C Mdiploid 0.00025 d 1 A 6	0.01199	-0.01355 – 0.03753	0.91995	0.35808612
C Mdiploid 0.00025 d 1 A 7	-0.00440	-0.03618 – 0.02738	-0.27152	0.78611237
C Mdiploid 0.00025 d 1 B 10	-0.00810	-0.04050 – 0.02431	-0.48987	0.62446295
C Mdiploid 0.00025 d 1 B 11	0.01723	-0.00711 – 0.04158	1.38737	0.16600725
C Mdiploid 0.00025 d 1 B 12	0.02984	-0.00604 – 0.06571	1.63013	0.10376368
C Mdiploid 0.00025 d 1 B 9	0.02155	-0.00464 – 0.04774	1.61298	0.10744003
C Mdiploid 0.00025 d 1 C 10	-0.01674	-0.05318 – 0.01970	-0.90042	0.36837253
C Mdiploid 0.00025 d 1 C 11	0.03217	0.00037 – 0.06396	1.98271	0.04799999
C Mdiploid 0.00025 d 1 C 12	0.01497	-0.01103 – 0.04097	1.12833	0.25977419
C Mdiploid 0.00025 d 1 C 9	-0.00518	-0.02846 – 0.01810	-0.43606	0.66299683
C Mdiploid 0.00025 d 1 D 2	0.00371	-0.10509 – 0.11251	0.06682	0.94675682
C Mdiploid 0.00025 d 1 E 2	0.01289	-0.00773 – 0.03352	1.22546	0.22103521
C Mdiploid 0.00025 d 1 H 2	0.00761	-0.02693 – 0.04215	0.43159	0.66624601
C Mdiploid 0.00025 d 1 H 3	-0.00459	-0.03880 – 0.02962	-0.26294	0.79271446
C Mdiploid 0.00025 d 1 H 4	0.01844	-0.02911 – 0.06599	0.75998	0.44765688
C Mdiploid 0.00025 d 1 H 5	0.02542	-0.00158 – 0.05242	1.84519	0.06565851
C Mdiploid 0.00025 d 1 H 6	0.00415	-0.02236 – 0.03065	0.30666	0.75924362
C Mdiploid 0.00025 d 1 H 7	0.00191	-0.05263 – 0.05644	0.06856	0.94537182
C Mdiploid 0.001 d 1 A 11	0.03612	0.00054 – 0.07170	1.98973	0.04721759
C Mdiploid 0.001 d 1 A 2	0.00809	-0.01969 – 0.03588	0.57084	0.56838732
C Mdiploid 0.001 d 1 A 2 B	0.00927	-0.01453 – 0.03307	0.76368	0.44545327
C Mdiploid 0.001 d 1 A 3	0.00032	-0.02879 – 0.02943	0.02166	0.98272673
C Mdiploid 0.001 d 1 A 3 B	0.00398	-0.02318 – 0.03114	0.28726	0.77404701
C Mdiploid 0.001 d 1 A 4	0.01781	-0.00788 – 0.04350	1.35859	0.17494851
C Mdiploid 0.001 d 1 A 4 B	0.00153	-0.02470 – 0.02777	0.11455	0.90885560
C Mdiploid 0.001 d 1 A 5	0.02938	-0.00138 – 0.06013	1.87221	0.06181742
C Mdiploid 0.001 d 1 A 5 B	-0.00512	-0.03385 – 0.02361	-0.34945	0.72690929
C Mdiploid 0.001 d 1 A 6	0.02327	-0.00226 – 0.04881	1.78643	0.07469382
C Mdiploid 0.001 d 1 A 6 B	0.00955	-0.01671 – 0.03581	0.71291	0.47626537
C Mdiploid 0.001 d 1 A 7	0.04291	0.01655 – 0.06926	3.19053	0.00151809
C Mdiploid 0.001 d 1 A 8	0.01513	-0.00780 – 0.03806	1.29286	0.19671558
C Mdiploid 0.001 d 1 A 9	0.02268	0.00057 – 0.04479	2.01065	0.04495102
C Mdiploid 0.001 d 1 B 1	0.05233	0.01476 – 0.08989	2.73013	0.00657598
C Mdiploid 0.001 d 1 B 2	0.04063	0.01391 – 0.06736	2.97975	0.00303894
C Mdiploid 0.001 d 1 C 1	0.05497	0.01942 – 0.09051	3.03106	0.00257568
C Mdiploid 0.001 d 1 C 2	0.03217	0.00671 – 0.05764	2.47643	0.01363235
C Mdiploid 0.001 d 1 D 4	0.08814	-0.01560 – 0.19188	1.66527	0.09654594
C Mdiploid 0.001 d 1 D 5	-0.01879	-0.04448 – 0.00690	-1.43345	0.15241540
C Mdiploid 0.001 d 1 D 6	0.04266	-0.04051 – 0.12583	1.00542	0.31522648
C Mdiploid 0.001 d 1 D 7	0.01221	-0.01136 – 0.03578	1.01531	0.31049575
C Mdiploid 0.001 d 1 D 8	-0.01271	-0.03390 – 0.00849	-1.17523	0.24051605
C Mdiploid 0.001 d 1 E 4	0.00302	-0.01233 – 0.01836	0.38541	0.70011148
C Mdiploid 0.001 d 1 E 5	0.03661	-0.01898 – 0.09221	1.29089	0.19739471
C Mdiploid 0.001 d 1 E 6	-0.02095	-0.03896 – -0.00293	-2.27831	0.02316980
C Mdiploid 0.001 d 1 E 7	0.02727	0.00695 – 0.04759	2.63081	0.00880654
C Mdiploid 0.001 d 1 E 8	-0.00435	-0.02578 – 0.01707	-0.39806	0.69077548
C Mdiploid 0.001 d 1 H 11	0.03103	0.00179 – 0.06027	2.07997	0.03808645
C Mdiploid 0.001 d 1 H 2	0.01207	-0.02316 – 0.04730	0.67164	0.50215015
C Mdiploid 0.001 d 1 H 2 B	-0.01772	-0.05147 – 0.01604	-1.02856	0.30423060
C Mdiploid 0.001 d 1 H 3	-0.00572	-0.04360 – 0.03215	-0.29622	0.76719449
C Mdiploid 0.001 d 1 H 3 B	0.00004	-0.03571 – 0.03580	0.00242	0.99807165
C Mdiploid 0.001 d 1 H 4	0.00887	-0.03993 – 0.05768	0.35635	0.72174082
C Mdiploid 0.001 d 1 H 4 B	-0.00269	-0.05249 – 0.04712	-0.10576	0.91582228
C Mdiploid 0.001 d 1 H 5	0.01359	-0.01637 – 0.04355	0.88908	0.37442760
C Mdiploid 0.001 d 1 H 5 B	-0.00275	-0.02940 – 0.02390	-0.20218	0.83986354
C Mdiploid 0.001 d 1 H 6	0.03459	0.00779 – 0.06138	2.53005	0.01174027
C Mdiploid 0.001 d 1 H 6 B	0.02557	-0.00048 – 0.05163	1.92371	0.05501270
C Mdiploid 0.001 d 1 H 7	0.00943	-0.04247 – 0.06132	0.35600	0.72200440
C Mdiploid 0.001 d 1 H 8	0.01372	-0.02957 – 0.05702	0.62116	0.53480274
C Mdiploid 0.001 d 1 H 9	0.01923	-0.02620 – 0.06465	0.82947	0.40727462
C Mdiploid 0.004 d 1 A 2	0.00074	-0.03253 – 0.03401	0.04350	0.96532350
C Mdiploid 0.004 d 1 A 3	0.01219	-0.02511 – 0.04950	0.64056	0.52212714
C Mdiploid 0.004 d 1 A 4	0.05056	0.01662 – 0.08451	2.91976	0.00367654
C Mdiploid 0.004 d 1 A 5	0.05149	0.01625 – 0.08674	2.86388	0.00437789
C Mdiploid 0.004 d 1 A 6	0.01932	-0.01445 – 0.05309	1.12129	0.26275545
C Mdiploid 0.004 d 1 A 7	0.02866	-0.01062 – 0.06794	1.42994	0.15341977
C Mdiploid 0.004 d 1 B 3	0.02463	-0.00370 – 0.05296	1.70377	0.08910621
C Mdiploid 0.004 d 1 B 4	0.03125	-0.00127 – 0.06378	1.88316	0.06031462
C Mdiploid 0.004 d 1 B 5	0.04299	0.01619 – 0.06979	3.14407	0.00177490
C Mdiploid 0.004 d 1 B 6	0.05328	0.01921 – 0.08734	3.06491	0.00230654
C Mdiploid 0.004 d 1 B 7	0.07750	0.04758 – 0.10742	5.07684	0.00000056
C Mdiploid 0.004 d 1 C 3	0.04274	0.01394 – 0.07154	2.90832	0.00381118
C Mdiploid 0.004 d 1 C 4	0.02312	-0.00410 – 0.05034	1.66497	0.09660494
C Mdiploid 0.004 d 1 C 5	0.04335	0.01800 – 0.06869	3.35220	0.00086846
C Mdiploid 0.004 d 1 C 6	0.06362	0.03622 – 0.09102	4.55122	0.00000685
C Mdiploid 0.004 d 1 C 7	0.03836	-0.01842 – 0.09514	1.32409	0.18613776
C Mdiploid 0.004 d 1 H 2	0.02225	-0.02299 – 0.06749	0.96402	0.33554719
C Mdiploid 0.004 d 1 H 3	0.02107	-0.02154 – 0.06367	0.96911	0.33300485
C Mdiploid 0.004 d 1 H 4	0.00548	-0.05419 – 0.06515	0.18007	0.85717391
C Mdiploid 0.004 d 1 H 5	0.06027	0.02281 – 0.09772	3.15346	0.00171997
C Mdiploid 0.004 d 1 H 6	0.04022	0.00560 – 0.07484	2.27725	0.02323376
C Mdiploid 0.004 d 1 H 7	0.01165	-0.05020 – 0.07350	0.36908	0.71224172
C Mhaploid 0.001 d 1 A 2	0.17748	-0.18514 – 0.54010	0.95929	0.33792087
C Mhaploid 0.001 d 1 A 3	0.03476	0.00257 – 0.06694	2.11667	0.03482853
C Mhaploid 0.001 d 1 A 3 B	0.02653	0.00424 – 0.04882	2.33283	0.02009125
C Mhaploid 0.001 d 1 A 4	0.07889	0.05620 – 0.10158	6.81600	<0.001
C Mhaploid 0.001 d 1 A 4 B	0.00801	-0.01550 – 0.03151	0.66754	0.50476447
C Mhaploid 0.001 d 1 A 5	0.07567	0.04879 – 0.10255	5.51778	0.00000006
C Mhaploid 0.001 d 1 A 5 B	0.05377	0.02908 – 0.07845	4.26875	0.00002393
C Mhaploid 0.001 d 1 A 6	0.02914	0.00850 – 0.04977	2.76747	0.00587909
C Mhaploid 0.001 d 1 A 6 B	0.02748	0.00967 – 0.04529	3.02424	0.00263330
C Mhaploid 0.001 d 1 A 7	0.28425	-0.16599 – 0.73448	1.23738	0.21658119
C Mhaploid 0.001 d 1 A 7 B	0.37549	-0.17001 – 0.92100	1.34912	0.17796797
C Mhaploid 0.001 d 1 H 2	0.05735	0.02810 – 0.08659	3.84317	0.00013869
C Mhaploid 0.001 d 1 H 3	0.00589	-0.02710 – 0.03888	0.34979	0.72665899
C Mhaploid 0.001 d 1 H 3 B	0.02272	-0.01348 – 0.05892	1.23012	0.21928455
C Mhaploid 0.001 d 1 H 4	0.06491	-0.01016 – 0.13999	1.69462	0.09083022
C Mhaploid 0.001 d 1 H 4 B	0.01455	-0.02446 – 0.05356	0.73088	0.46522961
C Mhaploid 0.001 d 1 H 5	0.24925	-0.16290 – 0.66141	1.18531	0.23651278
C Mhaploid 0.001 d 1 H 5 B	0.07769	0.04517 – 0.11021	4.68240	0.00000375
C Mhaploid 0.001 d 1 H 6	0.10010	0.04219 – 0.15802	3.38783	0.00076555
C Mhaploid 0.001 d 1 H 6 B	0.11257	0.06724 – 0.15790	4.86715	0.00000156
C Mhaploid 0.001 d 1 H 7	0.08866	0.04780 – 0.12952	4.25293	0.00002561
C Mhaploid 0.001 d 1 H 7 B	0.07939	0.04151 – 0.11728	4.10726	0.00004746
Observations	608
R2 / R2 adjusted	0.719 / 0.625

# -----------------------------------------------------------------------------

# create results dataframe
# ----------------------------------------------------
myres <- as.data.frame(coef(summary(mymod)), stringsAsFactors = F)
myres$p1s_l <- pt(coef(summary(mymod))[, 3], mymod$df, lower = T)  # one-sided test p-values for lower
myres$p1s_h <- pt(coef(summary(mymod))[, 3], mymod$df, lower = F)  # one-sided test p-values for higher
myres$p1s_l.adj <- p.adjust(myres$p1s_l, method = "fdr")  # fdr adjusted p value (for 1-sided lower)
myres$p1s_h.adj <- p.adjust(myres$p1s_h, method = "fdr")  # fdr adjusted p value (for 1-sided higher)
myres$bcid <- gsub("^.*?D", "", rownames(myres))  # fix barcode id names
myres$m_deltafit <- ddply(my, ~bctID, function(my) weighted.mean(my$deltafit, 
    my$rdeltafit))[, 2]
myres$m_rdeltafit <- ddply(my, ~bctID, function(my) mean(my$rdeltafit))[, 2]
myres$mcm_deltafit <- myres$m_deltafit - weighted.mean(myres$m_deltafit, myres$m_rdeltafit)
myres$color <- "fdr > 0.01"
myres$color[myres$p1s_h.adj <= 0.01] <- "fdr <= 0.01"
myres$color <- as.factor(myres$color)
myres$color <- relevel(myres$color, "fdr > 0.01")
myres$treatment <- gsub("0025d", "0025", gsub("001d", "001", gsub(".{4}$", "", 
    myres$bcid)))  # strip bc id from treatment label, store in myres
# -----------------------------------------------------------------------------

# review results
# --------------------------------------------------------------
paste0("weighted mean fitness increase for increasers = ", weighted.mean(myres$m_deltafit[myres$p1s_h.adj <= 
    0.01], myres$m_rdeltafit[myres$p1s_h.adj <= 0.01]))  # weighted mean fit inc of sig increasers only
## [1] "weighted mean fitness increase for increasers = 0.0679894955523455"
paste0("minimum fitness increase for increasers = ", min(myres$m_deltafit[myres$p1s_h.adj <= 
    0.01]))  # minimum significant increase of 2.4%
## [1] "minimum fitness increase for increasers = 0.0274814710527647"
paste0("maximum fitness increase for increasers = ", max(myres$m_deltafit[myres$p1s_h.adj <= 
    0.01]))  # maximum significant increase of 23.5%
## [1] "maximum fitness increase for increasers = 0.234510778858253"

# all treatments together
paste0("Number of barcodes with significant decreases in fitness = ", nrow(myres[myres$p1s_l.adj <= 
    0.01, ]))  # check for significant fitness decrease: # find none
## [1] "Number of barcodes with significant decreases in fitness = 0"
paste0("Number of barcodes with significant increases in fitness = ", nrow(myres[myres$p1s_h.adj <= 
    0.01, ]))  # check for significant fitness increase: # find 35 rows. (35 increasers out of 152)
## [1] "Number of barcodes with significant increases in fitness = 35"
paste0("Proportion of barcodes with significant increases in fitness = ", nrow(myres[myres$p1s_h.adj <= 
    0.01, ])/nrow(myres))  # prop increasers  # 0.23
## [1] "Proportion of barcodes with significant increases in fitness = 0.230263157894737"

# CM diploid 1:1000
paste0("CM diploid 1:1000:")
## [1] "CM diploid 1:1000:"
myrest <- myres[myres$treatment == "CMdiploid0.001", ]
paste0("n barcodes in treatment = ", nrow(myrest))
## [1] "n barcodes in treatment = 42"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 3"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.0714285714285714"

# CM haploid 1:1000
paste0("CM haploid 1:1000:")
## [1] "CM haploid 1:1000:"
myrest <- myres[myres$treatment == "CMhaploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 10"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.454545454545455"

# CM+salt diploid 1:1000
paste0("CM+salt diploid 1:1000:")
## [1] "CM+salt diploid 1:1000:"
myrest <- myres[myres$treatment == "CM+Saltdiploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 11"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.5"

# CM+ethanol diploid 1:1000
paste0("CM+ethanol diploid 1:1000:")
## [1] "CM+ethanol diploid 1:1000:"
myrest <- myres[myres$treatment == "CM+Ethanoldiploid0.001", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 2"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.0909090909090909"

# CM diploid 1:250
paste0("CM diploid 1:250:")
## [1] "CM diploid 1:250:"
myrest <- myres[myres$treatment == "CMdiploid0.004", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 9"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0.409090909090909"

# CM diploid 1:4000
paste0("CM diploid 1:4000:")
## [1] "CM diploid 1:4000:"
myrest <- myres[myres$treatment == "CMdiploid0.00025", ]
paste0("n entries in treatment = ", nrow(myrest))
## [1] "n entries in treatment = 22"
paste0("n barcodes with fitness increase in treatment = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ]))
## [1] "n barcodes with fitness increase in treatment = 0"
paste0("proportion of barcodes in treatment that increased in fitness = ", nrow(myrest[myrest$p1s_h.adj <= 
    0.01, ])/nrow(myrest))
## [1] "proportion of barcodes in treatment that increased in fitness = 0"
# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- myres$m_deltafit
myr <- myres$m_rdeltafit
myc <- myres$color
mybinwidth <- 0.01

mytitle <- NULL
myylab <- "Count"
myxlab <- "Fitness Change"

myhighlight1 <- "black"
myhighlight2 <- setpalette[1]
myhighlight3 <- setpalette[2]
myfillapha <- 0.25
myptalpha <- 0.35
mylinealpha <- 0.35
mypalette <- setpalette

mytextsize <- 8
myheight <- 3.3
mywidth <- 3.3


p <- ggplot(myres, aes(x = myx, color = myc, fill = myc)) + geom_histogram(alpha = myfillapha, 
    position = "identity", binwidth = mybinwidth) + geom_vline(xintercept = weighted.mean(myx, 
    myr), linetype = "dashed", color = myhighlight) + geom_vline(xintercept = weighted.mean(myx[myc == 
    "fdr > 0.01"], myr[myc == "fdr > 0.01"]), linetype = "dashed", color = myhighlight2) + 
    geom_vline(xintercept = weighted.mean(myx[myc == "fdr <= 0.01"], myr[myc == 
        "fdr <= 0.01"]), linetype = "dashed", color = myhighlight3) + geom_rug(alpha = myptalpha) + 
    scale_color_manual(values = mypalette, guide = FALSE) + scale_fill_manual(values = mypalette) + 
    theme_classic() + theme(legend.position = c(0.8, 0.8)) + labs(fill = "Significance") + 
    ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + geom_hline(yintercept = -0.01, 
    colour = "white", size = 1.01) + theme(axis.title.y = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = mytextsize, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = mytextsize, 
    angle = 0)) + theme(axis.text.x = element_text(size = mytextsize, angle = 0)) + 
    theme(legend.title = element_text(size = mytextsize, face = "bold")) + theme(legend.text = element_text(size = mytextsize))
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_3.pdf", width = mywidth, height = myheight)
p
dev.off()  # 7by7
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, myres, myrest, p)

Analysis: Assess power to detect fitness change in this experimental evolution system, given replicate variance in the POC fitness assay dataset.

# create a power table for lms
# -------------------------------------------------
myrmse <- 1.7568866199638  # this is the rmse from the fitness equivalence model*100 --> effect size of 1 corresonds to an ~1.76 fitness change 
myeffectsize <- 1/myrmse  # effect size of 1 is a 1.76 fitnesss change, so for a fitness change of 1.00, expect effect size of 1/myrmse
mypower.lm <- as.data.frame(cbind(rep(NA, times = 12), NA, NA, NA, NA))  # create dataframe
colnames(mypower.lm) <- c("u", "v", "f2", "sig.level", "power")  # fix column names
mypower.lm$u <- 1  # assign u
mypower.lm$v <- 6  # assign v
mypower.lm$sig.level <- c(0.05, 0.05, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01, 0.005, 
    0.005, 0.005, 0.005)  # assign significance
mypower.lm$f2 <- c(myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize * 
    5, myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize * 5, 
    myeffectsize * 0.1, myeffectsize, myeffectsize * 2, myeffectsize * 5)  # Assign effect size
for (i in 1:nrow(mypower.lm)) {
    # populate power column
    mytemp <- pwr.f2.test(u = mypower.lm$u[i], v = mypower.lm$v[i], sig.level = mypower.lm$sig.level[i], 
        f2 = mypower.lm$f2[i])
    mypower.lm$power[i] <- mytemp$power
}
# ------------------------------------------------------------------------------

# view and save
# ----------------------------------------------------------------
write.csv(mypower.lm, file = "000_Table_S5.csv")
# ------------------------------------------------------------------------------

rm(mytemp, i, mypower.lm)

Analysis: Assess effects of evoluitonary treatment on change in fitness in 250-generations of experimental evolution.

my <- myfa[order(myfa$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ----------------------------------------------------
mymod <- lmer(scale(deltafit) ~ scale(ccdeltafit) + treatment + (1 | bctID), 
    weights = rdeltafit, data = my, control = lmerControl(optimizer = "bobyqa", 
        optCtrl = list(maxfun = 20000)))
summary(mymod)  # problems with convergence here...solved with an optimizer and some extra iterations.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: scale(deltafit) ~ scale(ccdeltafit) + treatment + (1 | bctID)
##    Data: my
## Weights: rdeltafit
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
## 
## REML criterion at convergence: 789.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1056 -0.4647  0.0192  0.6301  3.5085 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  bctID    (Intercept)   0.1181  0.3437 
##  Residual             367.7128 19.1758 
## Number of obs: 608, groups:  bctID, 152
## 
## Fixed effects:
##                                  Estimate Std. Error        df t value
## (Intercept)                      -0.31415    0.05917 127.88336  -5.310
## scale(ccdeltafit)                -0.27053    0.01860 576.86686 -14.548
## treatmentCM+Ethanoldiploid0.001  -0.29452    0.10043 127.81225  -2.932
## treatmentCM+Saltdiploid0.001      0.96163    0.10787 145.79857   8.915
## treatmentCMdiploid0.00025        -0.08752    0.10119 129.15660  -0.865
## treatmentCMdiploid0.004           0.34621    0.10249 138.20327   3.378
## treatmentCMhaploid0.001           0.27966    0.10888 134.62106   2.568
##                                             Pr(>|t|)    
## (Intercept)                      0.00000047105690631 ***
## scale(ccdeltafit)               < 0.0000000000000002 ***
## treatmentCM+Ethanoldiploid0.001             0.003987 ** 
## treatmentCM+Saltdiploid0.001     0.00000000000000184 ***
## treatmentCMdiploid0.00025                   0.388701    
## treatmentCMdiploid0.004                     0.000949 ***
## treatmentCMhaploid0.001                     0.011306 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) scl(c) tCM+E0 tCM+S0 tCM0.000 tCM0.004
## scl(ccdltf) -0.020                                       
## trCM+E0.001 -0.592  0.164                                
## trCM+S0.001 -0.546 -0.091  0.308                         
## trCM0.00025 -0.585  0.016  0.347  0.319                  
## trtmCM0.004 -0.576 -0.060  0.330  0.322  0.336           
## trtmCM0.001 -0.545  0.112  0.338  0.288  0.319    0.307
# mymod <- lm(my$deltafit ~ my$ccdeltafit + my$treatment, weights =
# my$rdeltafit);r.squaredGLMM(mymod) # get the r2 for the lm version of the
# model if desired.

hist(residuals(mymod), breaks = 100)

hist(expand_residuals(mymod, my$rdeltafit), breaks = 100)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

anova(mymod)  # test term significance

# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------=
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    show.r2 = F, file = "000_Table_S8.html")  # r2 on the lm version of htis model is 0.3059, adjusted = 0.299

	scale(deltafit)
Predictors	Estimates	CI	Statistic	p
(Intercept)	-0.31415	-0.43012 – -0.19819	-5.30957	0.00000047
scale(ccdeltafit)	-0.27053	-0.30697 – -0.23408	-14.54794	<0.001
CM+Ethanoldiploid 0.001	-0.29452	-0.49136 – -0.09767	-2.93240	0.00398690
CM+Saltdiploid 0.001	0.96163	0.75021 – 1.17306	8.91470	<0.001
C Mdiploid 0.00025	-0.08752	-0.28584 – 0.11081	-0.86489	0.38870062
C Mdiploid 0.004	0.34621	0.14533 – 0.54710	3.37792	0.00094905
C Mhaploid 0.001	0.27966	0.06625 – 0.49307	2.56846	0.01130558
Random Effects
σ2	367.71
τ00 bctID	0.12
ICC	0.00
N bctID	152
Observations	608

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
mymod <- lmer((deltafit) ~ (ccdeltafit) + treatment + (1 | bctID), weights = rdeltafit, 
    data = my, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))  # problems with convergence here...solved with an optimizer and some extra iterations.
## Warning: Some predictor variables are on very different scales: consider
## rescaling

## Warning: Some predictor variables are on very different scales: consider
## rescaling
myx <- my$treatment
myy <- my$deltafit * 100
myr <- my$rdeltafit
myc <- my$treatment

mytitle <- NULL
myylab <- "Fitness Change"
myxlab <- "Treatment"

myoutline <- "grey20"
myfill <- "grey90"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(2, 3, 5, 6), y = c(14, 55, 15, 46), label = c("**", "***", "***", 
        "*"), size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_4.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, weights)

Analysis: Calculate fixation rate summaries by treatment for the 250-generation experiment.

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# Review Results
# -------------------------------------------------------------- my[my$fixed
# == T,] # these are the fixed wells my$treatment[my$fixed == T] # these are
# the treatments of the fixed wells

paste0("proportion of wells that attained fixation (n=76 total wells)= ", nrow(my[my$fixed == 
    T, ])/nrow(my))  # prop fixed
## [1] "proportion of wells that attained fixation (n=76 total wells)= 0.118421052631579"

# CM diploid 1:1000
myt <- my[my$treatment == "CMdiploid0.001", ]
paste0("n wells in CMdiploid0.001 = ", nrow(myt))  #  wells
## [1] "n wells in CMdiploid0.001 = 21"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 1"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0.0476190476190476"

# CM haploid 1:1000
myt <- my[my$treatment == "CMhaploid0.001", ]
paste0("n wells in CMhaploid0.001 = ", nrow(myt))  #  wells
## [1] "n wells in CMhaploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 5"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0.454545454545455"

# CM+Salt diploid 1:1000
myt <- my[my$treatment == "CM+Saltdiploid0.001", ]
paste0("n wells in CM+Saltdiploid0.001 = ", nrow(myt))  #  wells
## [1] "n wells in CM+Saltdiploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 3"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0.272727272727273"

# CM+Ethanol 1:1000
myt <- my[my$treatment == "CM+Ethanoldiploid0.001", ]
paste0("n wells in CM+Ethanoldiploid0.001 = ", nrow(myt))  #  wells
## [1] "n wells in CM+Ethanoldiploid0.001 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0"

# CM diploid 1:250
myt <- my[my$treatment == "CMdiploid0.004", ]
paste0("n wells in CMdiploid0.004 = ", nrow(myt))  #  wells
## [1] "n wells in CMdiploid0.004 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0"

# CM diploid 1:4000
myt <- my[my$treatment == "CMdiploid0.00025", ]
paste0("n wells in CMdiploid0.00025 = ", nrow(myt))  #  wells
## [1] "n wells in CMdiploid0.00025 = 11"
paste0("number of fixed wells in treatment = ", nrow(myt[myt$fixed, ]))  # number fixed
## [1] "number of fixed wells in treatment = 0"
paste0("proportion of fixed wells in treatment = ", nrow(myt[myt$fixed, ])/nrow(myt))  # prop fixed
## [1] "proportion of fixed wells in treatment = 0"
# -----------------------------------------------------------------------------
rm(my, myt)

Analysis: Assess effects of evolutionary treatment on T-MAX (generation of maximum deviation from generation-0 barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxdev ~ my$treatment
# + my$ccmaxdev_t0 + my$diff0, weights = my$rmaxdev_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$tmaxdev ~ my$treatment +
# my$diff0, weights = my$rmaxdev_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.

# without: mymod <- glm(my$tmaxdev ~ my$treatment + my$ccmaxdev, weights =
# my$rmaxdev); summary(mymod) # drop c as least significant term
mymod <- lm(my$tmaxdev ~ my$treatment, weights = my$rmaxdev)
summary(mymod)  # FINAL MODEL!
## 
## Call:
## lm(formula = my$tmaxdev ~ my$treatment, weights = my$rmaxdev)
## 
## Weighted Residuals:
##    Min     1Q Median     3Q    Max 
## -45513  -3616   1146   3910  16785 
## 
## Coefficients:
##                                    Estimate Std. Error t value
## (Intercept)                        233.5992    14.7420  15.846
## my$treatmentCM+Ethanoldiploid0.001 -20.8306    17.1490  -1.215
## my$treatmentCM+Saltdiploid0.001    -48.7017    18.3948  -2.648
## my$treatmentCMdiploid0.00025         2.2063    20.8325   0.106
## my$treatmentCMdiploid0.004         -43.3530    27.8363  -1.557
## my$treatmentCMhaploid0.001           0.8156    16.9681   0.048
##                                               Pr(>|t|)    
## (Intercept)                        <0.0000000000000002 ***
## my$treatmentCM+Ethanoldiploid0.001               0.229    
## my$treatmentCM+Saltdiploid0.001                  0.010 *  
## my$treatmentCMdiploid0.00025                     0.916    
## my$treatmentCMdiploid0.004                       0.124    
## my$treatmentCMhaploid0.001                       0.962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9517 on 70 degrees of freedom
## Multiple R-squared:  0.1972, Adjusted R-squared:  0.1398 
## F-statistic: 3.438 on 5 and 70 DF,  p-value: 0.007779
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rmaxdev), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.19716
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.1398143
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S9.html")

	my$tmaxdev
Predictors	Estimates	CI	Statistic	p
(Intercept)	233.59916	204.70539 – 262.49293	15.84584	<0.001
CM+Ethanoldiploid 0.001	-20.83060	-54.44209 – 12.78089	-1.21468	0.22856965
CM+Saltdiploid 0.001	-48.70167	-84.75481 – -12.64854	-2.64758	0.01000874
C Mdiploid 0.00025	2.20631	-38.62458 – 43.03720	0.10591	0.91595897
C Mdiploid 0.004	-43.35303	-97.91111 – 11.20506	-1.55743	0.12387850
C Mhaploid 0.001	0.81564	-32.44114 – 34.07241	0.04807	0.96179819
Observations	76
R2 / R2 adjusted	0.197 / 0.140

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxdev
myr <- my$rmaxdev
myc <- my$treatment

mytitle <- NULL
myylab <- "T-MAX (Generation)"
myxlab <- "Treatment"

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    # geom_hline(yintercept = 0, lty = 'dotted') +
coord_flip() + theme_classic() + theme(legend.position = "none") + ggtitle(mytitle) + 
    ylab(myylab) + xlab(myxlab) + scale_x_discrete(labels = c("CM \n Diploid \n 0.001", 
    "CM+EtOH \n Diploid \n 0.001", "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", 
    "CM \n Diploid \n 0.004", "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(3), y = c(255), label = "*", size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S4.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on M-MAX (magnitude of maximum deviation from generation-0 barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxdev ~ my$treatment
# + my$ccmaxdev_t0 + my$diff0, weights = my$rmaxdev_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$mmaxdev ~ my$treatment +
# my$diff0, weights = my$rmaxdev_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.

# without: mymod <- glm(my$mmaxdev ~ my$treatment + my$ccmaxdev, weights =
# my$rmaxdev); summary(mymod) # drop c as least significant term
mymod <- lm(my$mmaxdev ~ my$treatment, weights = my$rmaxdev)
summary(mymod)  # FINAL MODEL!
## 
## Call:
## lm(formula = my$mmaxdev ~ my$treatment, weights = my$rmaxdev)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -80.362 -11.841  -1.923  18.199 121.427 
## 
## Coefficients:
##                                    Estimate Std. Error t value  Pr(>|t|)
## (Intercept)                         0.22746    0.05437   4.184 0.0000820
## my$treatmentCM+Ethanoldiploid0.001 -0.06495    0.06325  -1.027   0.30799
## my$treatmentCM+Saltdiploid0.001     0.29240    0.06784   4.310 0.0000524
## my$treatmentCMdiploid0.00025        0.11117    0.07683   1.447   0.15240
## my$treatmentCMdiploid0.004          0.03904    0.10266   0.380   0.70491
## my$treatmentCMhaploid0.001          0.17791    0.06258   2.843   0.00585
##                                       
## (Intercept)                        ***
## my$treatmentCM+Ethanoldiploid0.001    
## my$treatmentCM+Saltdiploid0.001    ***
## my$treatmentCMdiploid0.00025          
## my$treatmentCMdiploid0.004            
## my$treatmentCMhaploid0.001         ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 35.1 on 70 degrees of freedom
## Multiple R-squared:  0.4564, Adjusted R-squared:  0.4176 
## F-statistic: 11.76 on 5 and 70 DF,  p-value: 0.00000002882
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rmaxdev), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.4564297
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.4176032
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S10.html")

	my$mmaxdev
Predictors	Estimates	CI	Statistic	p
(Intercept)	0.22746	0.12089 – 0.33402	4.18350	0.00008197
CM+Ethanoldiploid 0.001	-0.06495	-0.18891 – 0.05901	-1.02693	0.30798766
CM+Saltdiploid 0.001	0.29240	0.15943 – 0.42537	4.30998	0.00005236
C Mdiploid 0.00025	0.11117	-0.03942 – 0.26175	1.44685	0.15240115
C Mdiploid 0.004	0.03904	-0.16218 – 0.24025	0.38025	0.70491018
C Mhaploid 0.001	0.17791	0.05526 – 0.30057	2.84299	0.00585276
Observations	76
R2 / R2 adjusted	0.456 / 0.418

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxdev
myr <- my$rmaxdev
myc <- my$treatment

mytitle <- NULL
myylab <- "M-MAX"
myxlab <- "Treatment"

myoutline <- "grey20"
myfill <- "grey90"
myfillapha <- 0.35
myptalpha <- 0.35
mylinealpha <- 0.35

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(3, 6), y = c(0.91, 0.765), label = c("***", "**"), size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S5.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on T-MAX-RATE (generation of maximum rate of change in barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxrate ~ my$treatment
# + my$ccmaxrate_t0 + my$diff0, weights = my$rmaxrate_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$tmaxrate ~ my$treatment +
# my$diff0, weights = my$rmaxrate_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.

# without: mymod <- glm(my$tmaxrate ~ my$treatment + my$ccmaxrate, weights =
# my$rmaxrate); summary(mymod) # drop c as least significant term
mymod <- lm(my$tmaxrate ~ my$treatment, weights = my$rmaxrate)
summary(mymod)  # Final model!
## 
## Call:
## lm(formula = my$tmaxrate ~ my$treatment, weights = my$rmaxrate)
## 
## Weighted Residuals:
##    Min     1Q Median     3Q    Max 
## -36500  -2877    985   3554  25478 
## 
## Coefficients:
##                                    Estimate Std. Error t value
## (Intercept)                         235.251      7.709  30.518
## my$treatmentCM+Ethanoldiploid0.001   11.439     10.407   1.099
## my$treatmentCM+Saltdiploid0.001     -19.045     12.176  -1.564
## my$treatmentCMdiploid0.00025          9.910     17.741   0.559
## my$treatmentCMdiploid0.004            6.309     11.672   0.541
## my$treatmentCMhaploid0.001          -52.717     11.261  -4.681
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## my$treatmentCM+Ethanoldiploid0.001                0.275    
## my$treatmentCM+Saltdiploid0.001                   0.122    
## my$treatmentCMdiploid0.00025                      0.578    
## my$treatmentCMdiploid0.004                        0.591    
## my$treatmentCMhaploid0.001                    0.0000135 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9066 on 70 degrees of freedom
## Multiple R-squared:  0.3847, Adjusted R-squared:  0.3408 
## F-statistic: 8.755 on 5 and 70 DF,  p-value: 0.000001732
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rmaxrate), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.3847431
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.3407962
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S11.html")

	my$tmaxrate
Predictors	Estimates	CI	Statistic	p
(Intercept)	235.25147	220.14297 – 250.35998	30.51820	<0.001
CM+Ethanoldiploid 0.001	11.43864	-8.95850 – 31.83578	1.09914	0.27547276
CM+Saltdiploid 0.001	-19.04522	-42.90918 – 4.81874	-1.56420	0.12228104
C Mdiploid 0.00025	9.90963	-24.86202 – 44.68129	0.55857	0.57823532
C Mdiploid 0.004	6.30932	-16.56811 – 29.18674	0.54053	0.59054432
C Mhaploid 0.001	-52.71656	-74.78791 – -30.64520	-4.68130	0.00001354
Observations	76
R2 / R2 adjusted	0.385 / 0.341

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxrate
myr <- my$rmaxrate
myc <- my$treatment

mytitle <- NULL
myylab <- "T-MAX-RATE"
myxlab <- "Treatment"

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    # geom_hline(yintercept = 0, lty = 'dotted') +
coord_flip() + theme_classic() + theme(legend.position = "none") + ggtitle(mytitle) + 
    ylab(myylab) + xlab(myxlab) + scale_x_discrete(labels = c("CM \n Diploid \n 0.001", 
    "CM+EtOH \n Diploid \n 0.001", "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", 
    "CM \n Diploid \n 0.004", "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(6), y = c(257), label = "***", size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S6.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on M-MAX-RATE (magnitude of maximum rate of change in barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxrate ~ my$treatment
# + my$ccmaxrate_t0 + my$diff0, weights = my$rmaxrate_rt0); summary(mymod) #
# drop c as least significant term mymod <- glm(my$mmaxrate ~ my$treatment +
# my$diff0, weights = my$rmaxrate_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.

# without: mymod <- glm(my$mmaxrate ~ my$treatment + my$ccmaxrate, weights =
# my$rmaxrate); summary(mymod) # drop c as least significant term
mymod <- lm(my$mmaxrate ~ my$treatment, weights = my$rmaxrate)
summary(mymod)  # FINAL MODEL!
## 
## Call:
## lm(formula = my$mmaxrate ~ my$treatment, weights = my$rmaxrate)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0806 -0.8103 -0.1020  0.8827  4.2056 
## 
## Coefficients:
##                                     Estimate Std. Error t value
## (Intercept)                         0.009582   0.001258   7.620
## my$treatmentCM+Ethanoldiploid0.001 -0.002035   0.001698  -1.199
## my$treatmentCM+Saltdiploid0.001     0.003154   0.001986   1.588
## my$treatmentCMdiploid0.00025        0.006305   0.002894   2.178
## my$treatmentCMdiploid0.004         -0.002101   0.001904  -1.104
## my$treatmentCMhaploid0.001         -0.002801   0.001837  -1.525
##                                           Pr(>|t|)    
## (Intercept)                        0.0000000000916 ***
## my$treatmentCM+Ethanoldiploid0.001          0.2347    
## my$treatmentCM+Saltdiploid0.001             0.1168    
## my$treatmentCMdiploid0.00025                0.0328 *  
## my$treatmentCMdiploid0.004                  0.2736    
## my$treatmentCMhaploid0.001                  0.1319    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.479 on 70 degrees of freedom
## Multiple R-squared:  0.2098, Adjusted R-squared:  0.1534 
## F-statistic: 3.718 on 5 and 70 DF,  p-value: 0.004837
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rmaxrate), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.2098238
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.1533827
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S12.html")

	my$mmaxrate
Predictors	Estimates	CI	Statistic	p
(Intercept)	0.00958	0.00712 – 0.01205	7.61961	<0.001
CM+Ethanoldiploid 0.001	-0.00204	-0.00536 – 0.00129	-1.19871	0.23468475
CM+Saltdiploid 0.001	0.00315	-0.00074 – 0.00705	1.58776	0.11684830
C Mdiploid 0.00025	0.00630	0.00063 – 0.01198	2.17830	0.03275392
C Mdiploid 0.004	-0.00210	-0.00583 – 0.00163	-1.10353	0.27357441
C Mhaploid 0.001	-0.00280	-0.00640 – 0.00080	-1.52455	0.13187664
Observations	76
R2 / R2 adjusted	0.210 / 0.153

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxrate
myr <- my$rmaxrate
myc <- my$treatment

mytitle <- NULL
myylab <- "M-MAX-RATE"
myxlab <- "Treatment"


mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(4), y = c(0.0325), label = "*", size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S7.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on T-MAX-DIFF (generation of maximum difference in sympatric barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tmaxdiff ~ my$treatment
# + my$ccmaxdiff_t0 + my$diff0, weights = my$rmaxdiff_rt0); summary(mymod) #
# drop c as least significant term
mymod <- lm(my$tmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
summary(mymod)  # # FINAL MODEL!
## 
## Call:
## lm(formula = my$tmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
## 
## Weighted Residuals:
##    Min     1Q Median     3Q    Max 
## -58426  -5302    574   8332  59655 
## 
## Coefficients:
##                                    Estimate Std. Error t value
## (Intercept)                         241.970     34.101   7.096
## my$treatmentCM+Ethanoldiploid0.001  -22.732     35.785  -0.635
## my$treatmentCM+Saltdiploid0.001      -5.533     38.822  -0.143
## my$treatmentCMdiploid0.00025         11.312     77.413   0.146
## my$treatmentCMdiploid0.004          -23.716     63.676  -0.372
## my$treatmentCMhaploid0.001            1.612     35.479   0.045
## my$diff0                           -328.617     81.224  -4.046
##                                          Pr(>|t|)    
## (Intercept)                        0.000000000892 ***
## my$treatmentCM+Ethanoldiploid0.001       0.527381    
## my$treatmentCM+Saltdiploid0.001          0.887074    
## my$treatmentCMdiploid0.00025             0.884245    
## my$treatmentCMdiploid0.004               0.710705    
## my$treatmentCMhaploid0.001               0.963902    
## my$diff0                                 0.000134 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19160 on 69 degrees of freedom
## Multiple R-squared:  0.2342, Adjusted R-squared:  0.1676 
## F-statistic: 3.517 on 6 and 69 DF,  p-value: 0.004283

anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rtcc), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.2342223
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.1676329
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S13.html")

	my\(tmaxdiff</th> </tr> <tr> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; text-align:left; ">Predictors</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Estimates</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">CI</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Statistic</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">p</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">241.97046</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">175.13346&nbsp;&ndash;&nbsp;308.80745</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">7.09567</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>&lt;0.001</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Ethanoldiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-22.73154</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-92.86866&nbsp;&ndash;&nbsp;47.40558</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.63523</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.52738134</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Saltdiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-5.53349</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-81.62393&nbsp;&ndash;&nbsp;70.55695</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.14253</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.88707398</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.00025</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">11.31234</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-140.41438&nbsp;&ndash;&nbsp;163.03905</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.14613</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.88424511</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.004</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-23.71583</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-148.51929&nbsp;&ndash;&nbsp;101.08763</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.37244</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.71070458</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mhaploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">1.61153</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-67.92574&nbsp;&ndash;&nbsp;71.14880</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.04542</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.96390191</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">my\)diff 0				-328.61667	-487.81355 – -169.41979	-4.04579	0.00013406
Observations	76
R2 / R2 adjusted	0.234 / 0.168

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tmaxdiff
myr <- my$rmaxdiff_rt0
myc <- my$treatment

mytitle <- NULL
myylab <- "T-MAX-DIFF"
myxlab <- "Treatment"

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0))
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S8.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on M-MAX-DIFF (magnitude of maximum difference in sympatric barcoded strain abundance)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$mmaxdiff ~ my$treatment
# + my$ccmaxdiff_t0 + my$diff0, weights = my$rmaxdiff_rt0); summary(mymod) #
# drop c as least significant term
mymod <- lm(my$mmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
summary(mymod)  # # FINAL MODEL!
## 
## Call:
## lm(formula = my$mmaxdiff ~ my$treatment + my$diff0, weights = my$rmaxdiff_rt0)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -90.008  -8.179   2.215  10.833  47.792 
## 
## Coefficients:
##                                    Estimate Std. Error t value   Pr(>|t|)
## (Intercept)                         0.18985    0.04004   4.741 0.00001105
## my$treatmentCM+Ethanoldiploid0.001  0.00454    0.04202   0.108   0.914271
## my$treatmentCM+Saltdiploid0.001     0.17472    0.04559   3.833   0.000277
## my$treatmentCMdiploid0.00025        0.15524    0.09090   1.708   0.092165
## my$treatmentCMdiploid0.004          0.03311    0.07477   0.443   0.659311
## my$treatmentCMhaploid0.001          0.20429    0.04166   4.904 0.00000601
## my$diff0                            0.29644    0.09538   3.108   0.002734
##                                       
## (Intercept)                        ***
## my$treatmentCM+Ethanoldiploid0.001    
## my$treatmentCM+Saltdiploid0.001    ***
## my$treatmentCMdiploid0.00025       .  
## my$treatmentCMdiploid0.004            
## my$treatmentCMhaploid0.001         ***
## my$diff0                           ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.49 on 69 degrees of freedom
## Multiple R-squared:  0.4789, Adjusted R-squared:  0.4336 
## F-statistic: 10.57 on 6 and 69 DF,  p-value: 0.00000002701
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rtcc), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.4789295
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.433619
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S14.html")

	my\(mmaxdiff</th> </tr> <tr> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; text-align:left; ">Predictors</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Estimates</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">CI</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">Statistic</td> <td style=" text-align:center; border-bottom:1px solid; font-style:italic; font-weight:normal; ">p</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">(Intercept)</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.18985</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.11137&nbsp;&ndash;&nbsp;0.26833</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">4.74128</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00001105</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Ethanoldiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.00454</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.07782&nbsp;&ndash;&nbsp;0.08690</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.10805</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.91427138</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">CM+Saltdiploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.17472</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.08537&nbsp;&ndash;&nbsp;0.26406</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">3.83269</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00027655</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.00025</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.15524</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.02292&nbsp;&ndash;&nbsp;0.33340</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">1.70782</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.09216521</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mdiploid 0.004</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.03311</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">-0.11344&nbsp;&ndash;&nbsp;0.17965</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.44278</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.65931134</td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">C Mhaploid 0.001</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.20429</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">0.12264&nbsp;&ndash;&nbsp;0.28594</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; ">4.90375</td> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:center; "><strong>0.00000601</strong></td> </tr> <tr> <td style=" padding:0.2cm; text-align:left; vertical-align:top; text-align:left; ">my\)diff 0				0.29644	0.10951 – 0.48337	3.10815	0.00273414
Observations	76
R2 / R2 adjusted	0.479 / 0.434

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$mmaxdiff
myr <- my$rmaxdiff_rt0
myc <- my$treatment

mytitle <- NULL
myylab <- "M-MAX-DIFF"
myxlab <- "Treatment"

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(3, 6), y = c(0.525), label = "***", size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S9.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)

Analysis: Assess effects of evolutionary treatment on total change in barcode abundance (summed across six intervals)

my <- myevo[order(myevo$bctID), ]  # placeholder dataframe
# model & diagnostic plots
# ---------------------------------------------------- with initial
# proportions included in the model: mymod <- glm(my$tcc ~ my$treatment +
# my$cctcc_t0 + my$diff0, weights = my$rtcc_rt0); summary(mymod) # drop c as
# least significant term mymod <- glm(my$tcc ~ my$treatment + my$diff0,
# weights = my$rtcc_rt0); summary(mymod) # diff 0 still not
# significant...run the model set without it.

# without: mymod <- glm(my$tcc ~ my$treatment + my$cctcc, weights =
# my$rtcc); summary(mymod) # drop c as least significant term
mymod <- lm(my$tcc ~ my$treatment, weights = my$rtcc)
summary(mymod)  # FINAL MODEL!
## 
## Call:
## lm(formula = my$tcc ~ my$treatment, weights = my$rtcc)
## 
## Weighted Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.7198 -1.3007 -0.0058  1.0657  5.4905 
## 
## Coefficients:
##                                      Estimate Std. Error t value
## (Intercept)                         0.0219194  0.0027409   7.997
## my$treatmentCM+Ethanoldiploid0.001 -0.0072401  0.0032633  -2.219
## my$treatmentCM+Saltdiploid0.001     0.0046752  0.0035042   1.334
## my$treatmentCMdiploid0.00025        0.0182730  0.0056216   3.251
## my$treatmentCMdiploid0.004         -0.0007325  0.0043115  -0.170
## my$treatmentCMhaploid0.001         -0.0066256  0.0033076  -2.003
##                                           Pr(>|t|)    
## (Intercept)                        0.0000000000185 ***
## my$treatmentCM+Ethanoldiploid0.001         0.02976 *  
## my$treatmentCM+Saltdiploid0.001            0.18647    
## my$treatmentCMdiploid0.00025               0.00177 ** 
## my$treatmentCMdiploid0.004                 0.86559    
## my$treatmentCMhaploid0.001                 0.04904 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.15 on 70 degrees of freedom
## Multiple R-squared:  0.3746, Adjusted R-squared:   0.33 
## F-statistic: 8.387 on 5 and 70 DF,  p-value: 0.000002951
anova(mymod)  # test term significance with an lm 

hist(residuals(mymod), breaks = 25)

hist(expand_residuals(mymod, my$rtcc), breaks = 25)  # this more accurately represents the distribution of residuals in the model output. -- tails look fine here. 

plot(mymod)  # check model fit 

summary(mymod)$r.squared  # for R2 value
## [1] 0.3746437
summary(mymod)$adj.r.squared  # for adjusted R2 value
## [1] 0.3299754
# -----------------------------------------------------------------------------

# model output
# -----------------------------------------------------------------
tab_model(mymod, digits = 5, digits.p = 8, show.stat = T, wrap.labels = 100, 
    file = "000_Table_S15.html")

	my$tcc
Predictors	Estimates	CI	Statistic	p
(Intercept)	0.02192	0.01655 – 0.02729	7.99719	<0.001
CM+Ethanoldiploid 0.001	-0.00724	-0.01364 – -0.00084	-2.21861	0.02975813
CM+Saltdiploid 0.001	0.00468	-0.00219 – 0.01154	1.33418	0.18646947
C Mdiploid 0.00025	0.01827	0.00725 – 0.02929	3.25053	0.00177244
C Mdiploid 0.004	-0.00073	-0.00918 – 0.00772	-0.16989	0.86559011
C Mhaploid 0.001	-0.00663	-0.01311 – -0.00014	-2.00311	0.04903884
Observations	76
R2 / R2 adjusted	0.375 / 0.330

# -----------------------------------------------------------------------------

# Visualize
# -------------------------------------------------------------------
myx <- my$treatment
myy <- my$tcc
myr <- my$rtcc
myc <- my$treatment

mytitle <- NULL
myylab <- "Total Change in Abundance"
myxlab <- "Treatment"

mypalette <- setpalette
mytextsize <- 8
myheight <- 3.5
mywidth <- 7
myptsize <- ((((myr/max(myr)) * 10) * ((myheight * mywidth)/(7^2))) - min(((myr/max(myr)) * 
    10) * ((myheight * mywidth)/(7^2)))) + 1


p <- ggplot(my, aes(x = myx, y = myy, weight = myr, color = myc)) + geom_violin(color = myoutline, 
    fill = myfill, draw_quantiles = NULL, inherit.aes = F, aes(x = myx, y = myy), 
    trim = T) + geom_jitter(width = 0.15, height = 0, alpha = myptalpha, size = myptsize) + 
    geom_segment(x = 5.9, xend = 6.1, y = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), yend = weighted.mean(myy[myx == "CMhaploid0.001"], 
        myr[myx == "CMhaploid0.001"]), color = myoutline) + geom_segment(x = 4.9, 
    xend = 5.1, y = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    yend = weighted.mean(myy[myx == "CMdiploid0.004"], myr[myx == "CMdiploid0.004"]), 
    color = myoutline) + geom_segment(x = 3.9, xend = 4.1, y = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), yend = weighted.mean(myy[myx == 
    "CMdiploid0.00025"], myr[myx == "CMdiploid0.00025"]), color = myoutline) + 
    geom_segment(x = 2.9, xend = 3.1, y = weighted.mean(myy[myx == "CM+Saltdiploid0.001"], 
        myr[myx == "CM+Saltdiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Saltdiploid0.001"], myr[myx == "CM+Saltdiploid0.001"]), color = myoutline) + 
    geom_segment(x = 1.9, xend = 2.1, y = weighted.mean(myy[myx == "CM+Ethanoldiploid0.001"], 
        myr[myx == "CM+Ethanoldiploid0.001"]), yend = weighted.mean(myy[myx == 
        "CM+Ethanoldiploid0.001"], myr[myx == "CM+Ethanoldiploid0.001"]), color = myoutline) + 
    geom_segment(x = 0.9, xend = 1.1, y = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), yend = weighted.mean(myy[myx == "CMdiploid0.001"], 
        myr[myx == "CMdiploid0.001"]), color = myoutline) + geom_rug(sides = "b", 
    position = position_nudge(x = -0.75), alpha = myptalpha) + scale_color_manual(values = mypalette) + 
    geom_hline(yintercept = 0, lty = "dotted") + coord_flip() + theme_classic() + 
    theme(legend.position = "none") + ggtitle(mytitle) + ylab(myylab) + xlab(myxlab) + 
    scale_x_discrete(labels = c("CM \n Diploid \n 0.001", "CM+EtOH \n Diploid \n 0.001", 
        "CM+NaCl \n Diploid \n 0.001", "CM \n Diploid \n 0.00025", "CM \n Diploid \n 0.004", 
        "CM \n Haploid \n 0.001")) + theme(axis.title.y = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.title.x = element_text(size = 8, 
    face = "bold", margin = margin(0, 0, 0, 0, "pt"))) + theme(axis.text.y = element_text(size = 8, 
    angle = 0)) + theme(axis.text.x = element_text(size = 8, angle = 0)) + annotate("text", 
    x = c(2, 4, 6), y = c(0.0275, 0.06825, 0.026), label = c("*", "**", "*"), 
    size = 8, fontface = "bold")
p

# -----------------------------------------------------------------------------

# figure output
# ---------------------------------------------------------------
pdf(file = "000_Figure_S10.pdf", width = mywidth, height = myheight)
p
dev.off()
## quartz_off_screen 
##                 2
# -----------------------------------------------------------------------------
rm(my, mymod, p)